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Figure  4:  Variation  of  Hall  probe  field  measurement  with  the  probe  on  a  uniform  (i.e. 
unflawed)  copper  plate  as  a  function  of  frequency,  comparison  of  measurements,  (circles) 
with  theory,  (solid  line). 


axial  sensor  were  carried  out  at  2500  Hz.  The  predicted  amplitude  and  phase  of  the  signal 
as  a  function  of  position  for  probe  P3  (Table  1)  is  shown  in  Figure  5  together  with  the 
measurements.  Note  that  the  noise  in  the  signal  is  hardly  discernable  on  the  graphs. 

The  predictions  are  somewhat  dependent  on  the  number  of  volume  cells  used,  although, 
with  a  sufficient  number  the  results  is  stable.  With  a  12x12x1  array  of  volume  elements 
the  simulation  code,  taking  5  minutes  to  run  on  a  1.5GHz  PC,  is  in  good  agreement  with 
experiment  both  for  the  amplitude  and  the  phase  data.  The  over-prediction  error  in  the 
amplitude  has  been  reduced  by  improving  the  calculation  and  the  calibration  precision. 
However,  in  view  of  the  accuracy  of  the  prediction  on  an  unflawed  plate,  which  is  within 
4%,  we  aim  in  the  future  to  reduce  the  discrepancy  for  flaws  to  around  5%. 

The  phase  discrepancy  for  the  recess  calculation  is  about  10  degrees.  This  is  not 
unreasonable  in  view  of  the  fact  the  phase  can  change  markedly  with  the  depth  of  the 
subsurface  flaw  or  with  frequency,  as  shown  in  Figure  6.  However  an  effort  will  be  made 
in  the  future  to  reduce  the  phase  discrepancy  further. 

1.2.3  Holes  and  fasteners 

One  of  the  aims  of  the  current  project  is  to  predict  signals  due  to  structural  features  such 
as  holes,  fasteners  and  lap  joints  and  then  perform  calculations  involving  these  features  in 
the  presence  of  material  loss  due  to  corrosion.  Work  on  the  combined  problem  has  not  yet 
been  completed  but  experiments  on  holes,  fasteners  and  lap  joints  have  been  performed. 

Experimental  measurements  of  the  magnetic  field  due  to  circular  holes  in  an  aluminum 
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Figure  6:  Theoretical  predictions  of  the  variation  of  magnetic  field  measurements  due  to  a 
square  recess  measurements  at  a  frequencies  0.5kHz  (data  1)  1.5kHz  (data  2)  and  2.5  kHz 
(data  3).  Note  that  the  phase  of  the  response  changes  by  over  100  degrees  for  a  frequency 
change  from  0.5kHz  to  2.5kHz. 

1.3  Conclusions 

•  Progress:  It  was  our  intention,  at  the  start  of  the  current  project,  to  validate  fully 
a  new  volume  element  code  based  on  a  discrete  representation  of  the  electromagnetic 
field  in  terms  of  edge  elements.  This  as  not  been  possible  firstly  because  the  original 
plan  was  perhaps  a  little  ambitious.  In  our  optimism  it  was  felt  that  a  new  volume 
element  formulation  could  be  developed  and  coded  within  a  few  months.  Other 
factors  included  the  delays  in  finalizing  the  contractual  arrangements  between  ISU 
and  Cassino  which  lead  to  delays  in  completion  of  the  formulation  thus  limiting  the 
development  time  that  could  be  spent  on  the  new  code.  At  this  stage,  the  edge 
element  formulation  has  been  coded  for  a  half-space  conductor  rather  than  a  layered 
conductor,  therefore  it  need  further  work  before  calculations  can  be  compared  with 
measurements  on  flaws  in  plates.  However  selected  comparisons  between  theory  and 
experiment  have  been  carried  out  using  the  earlier  code  based  on  a  traditional  volume 
element  scheme. 

•  Sensors:  Experiments  have  been  performed  with  a  single  sensor  probe  to  test  the¬ 
oretical  results  and  to  refine  experimental  techniques.  Two  precision  made  probes 
with  Hall  sensors  have  been  built  to  perform  the  measurements  although  in  fact,  all 
the  experimental  data  presented  here  is  for  a  probe  with  a  Honeywell  634SS2  sen¬ 
sor.  The  Honeywell  was  first  choice  simply  because  it  is  familiar  to  us.  However 
the  Asahi  Kasia  HW-108A  has  a  good  specification  on  paper,  has  the  advantage  of 
being  enclosed  in  a  small  package  and  therefore  seems  suitable  for  a  high  density 
array.  A  single  line  of  such  sensors  could  be  mounted  less  than  2.5  mm  apart.  By 
using  say  three  parallel  lines  of  sensor  in  a  staggered  array,  it  would  be  possible  to 
get  measurements  in  the  direction  of  the  senor  lines  with  a  separation  of  about  than 
1  mm. 
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Figure  7:  Magnetic  field  measurements  of  holes  in  a  conducting  plate. 


•  Comparisons:  The  experimental  technique  has  been  greatly  improved  during  this 
phase  of  the  project  and  there  is  no  doubt  that  the  data  is  calibrated  correctly  and  is 
reliable.  Overall  the  comparison  between  theory  and  experiment  are  reasonable  and 
usually  within  10%.  The  discrepancies  are  thought  to  be  related  to  limitations  of  the 
model  rather  than  errors  in  the  experiment.  It  is  to  be  hope  that  the  edge  element 
model  will  provide  further  improvements  in  the  accuracy  of  the  predictions. 

•  Imaging:  While  awaiting  the  completion  of  the  new  formulation  by  colleagues  at 
the  University  of  Cassino,  the  question  of  producing  quantitative  images  of  a  corrode 
region  was  address.  This  topic  is  taken  up  in  the  following  article.  At  the  time  of 
writing,  no  numerical  experiment  have  not  been  perform  to  illustrate  the  use  of  the 
imaging  technique  but  given  the  simplicity  of  the  approach,  this  will  not  be  difficult 
and  should  be  carried  out  in  the  near  future. 
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Figure  11:  Magnetic  field  predictions  for  a  simplified  lap  joint  45  mm  wide  forming  a  2 
plate  stack  with  a  0.14  mm  gap  between  plates.  The  in-phase  (real  part)  and  quadrature 
component  (imaginary  part)  of  the  field  at  the  axial  sensor  are  shown  as  a  function  of 
position.  The  upper  pair  of  results  are  for  a  frequency  of  500Hz  and  the  lower  pair  for  a 
frequencies  of  2.5  kHz. 


Magnetic  Field  at  a  Lap  Joint  Magnetic  Field  at  a  Lap  Joint 


Figure  12:  Comparison  of  the  predicted  and  measured  axial  magnetic  field  variation  with 
probe  position  above  a  lap  joint.  The  field  component  in  phase  with  the  coil  current  and 
in  quadrature  with  the  coil  current  are  shown.  These  are,  the  real  and  imaginary  parts 
respectively  of  the  phasor  representing  the  magnetic  field. 
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F2=  £  \Htpt)  -  H(cm)l2  (A2) 

m=l,M 

The  Hall  device  liftoff  value,  zh,  is  them  varied  until  this  error  function  reaches  a 
minimum. 


FIGURE  Al.  Probe  cross  section  showing  key  probe  parameters. 
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by  the  driver.  A  plate  or  a  system  of  parallel  plates  perpendicular  to  the  z-axis  has 
a  translational  invariance  which  allowing  the  unperturbed  electric  field  to  be  written  as 
E(«)(x  —  x',y  —  y',z')  with  x  and  y  as  the  probe  coordinates.  For  pick-up  and  driver, 
n  =  1,2  respectively. 

For  small  amount  of  material  loss  due  to  corrosion  at  the  surface  of  a  plate,  the  inte¬ 
gration  over  the  direction  normal  to  the  surface  of  the  unflawed  plate,  the  z— direction  in 
this  case,  can  be  approximated  as 

J  E(0)(r')  •  P(r ')dz'  =  cr0E(1)(a;  -  x',y  -  y',z0)  ■  E(2)(x  -  x',y-  y' ,  z0)5z(x' ,y')  (2.6) 

where  5z(x',yr )  is  the  height  of  the  loss  region,  zq  is  the  z— coordinate  of  the  the  original 
(unflawed)  surface  and  E^(x7,  y' ,  zq)  is  the  unperturbed  field  due  to  the  driver  coil.  This 
approximation  can  be  justified  by  using  the  approach  described  in  reference  [1]  where  an 
expression  for  the  variation  of  the  probe  response  is  given  for  a  small  variation  in  the  loca¬ 
tion  of  a  conducting  surface.  The  variation  Sz(x',  y')  is  interpreted  as  being  associated  with 
material  loss  due  to  corrosion  and  the  variation  of  probe  signal,  written  as  V12,  represents 
the  signal  due  to  material  loss. 

Substituting  (2.6)  into  (2.5)  and  putting 

HV2{x  -x',y-  y)  =  -o-0E(1)(z  -x',y-  y  ,  z0)  •  E{2)(x  -x',y-  y  ,  z0)  (2.7) 

gives 

/oo  poo 

/  H(x-x',y-y')8z(x' ,y')dx' dy'  (2.8) 

-00  J  —  OO 

which  is  a  two  dimensional  convolution  for  the  material  loss  8z(x,y). 

2.4  Conclusion 

The  thesis  presented  here  is  that  the  deconvolution  via  Weiner  filtering  applied  to  eddy 
current  data  images  of  surface  cracks  has  not  been  justified  theoretically  and  is  not  effective. 
However,  the  same  procedure  can  be  used  to  find  the  material  loss  from  a  flat  plate  in  a 
stack  of  flat  plates.  The  basic  requirement  is  that  the  loss  estimator  5z(x,  y),  must  be  small 
compared  to  the  plate  thickness. 

A  key  practical  question  for  further  consideration  is  one  of  resolution.  It  is  well  known 
that  corrosion  can  be  intensely  localize  in  the  form  of  pits  which  may  be  small.  The 
ability  to  find  the  depth  of  such  features  using  the  Fast  Fourier  Transform  methods  may  be 
limited  because  of  both  the  fundamental  limitation  of  the  eddy  current  modality,  it  is  after 
all  a  highly  diffusive  phenomena,  and  in  addition  there  are  sampling  limits  which  reduce 
resolution.  The  former  can  be  mitigated  by  using  higher  frequencies,  the  latter  by  using  a 
high  density  sensor  array. 
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Article  3 

Eddy  current  excitation  using  a  racetrack  coil  with 
a  sensor  array  for  magnetic  field  measurement 

J.  R.  Bowler  and  V.  Katyal 

Iowa  State  University,  Center  Nondestructive  Evaluation, 

ASC  II,  1915  Scholl  Road,  Ames  Iowa  50011. 

Abstract  Calculations  have  been  performed  to  determine  the  response 
of  a  new  eddy  current  probe  for  the  detection  of  subsurface  flaws  in 
planar  multilayered  structures.  The  probe  consists  of  a  racetrack  coil 
and  a  linear  array  of  solid  state  sensors  for  detecting  perturbations  in 
the  electromagnetic  field  due  to  defects.  The  sensor  array  allows  field 
measurements  to  be  made  at  a  number  of  closely  spaced  locations  with¬ 
out  moving  the  probe  and  thereby  accelerates  the  inspection  process.  A 
magnetic  shell  model  of  the  probe  is  used  for  finding  the  electric  field  in 
the  unflawed  structure.  The  fields  due  to  the  linear  “straights”  and  the 
semicircular  “bends”  are  found  separately  and  added  to  give  the  com¬ 
bined  field  of  the  racetrack  coil.  The  flaw  response  is  then  computed 
using  a  volume  element  calculation.  In  order  to  validate  the  calculation, 
field  predictions  for  a  racetrack  coil  having  straights  of  zero  length  are 
compared  with  results  for  a  circular  coil.  The  results  are  found  to  be 
consistent. 

3.1  Introduction 

In  eddy  current  inspection,  an  induction  coil  is  often  used  both  to  induce  current  in  a 
conducting  component  and  to  detect  magnetic  field  perturbations  due  to  flaws.  For  sub¬ 
surface  defects,  a  low  frequency  excitation  ensures  an  adequate  depth  of  field  penetration. 
However,  at  lower  frequencies  the  effectiveness  of  the  coil  as  both  inducer  and  sensor  is 
diminished  since  electromagnetic  induction  depends  on  the  rate  of  change  of  magnetic  flux. 
To  overcome  the  limitations  of  the  induction  coil  as  a  low  frequency  field  sensor,  a  solid 
state  device,  such  as  a  giant  magneto-resistor  or  Hall  sensor,  can  be  used  instead.  A  coil 
used  only  as  driver  can  be  larger  than  otherwise  without  compromising  the  spatial  reso¬ 
lution  of  the  measurements.  The  large  coil  can  produce  a  greater  field  while  good  spatial 
resolution  is  obtained  by  using  small  sensors. 

This  article  gives  the  analysis  of  an  eddy  current  probe  for  the  detection  of  subsurface 
flaws  in  multilayered  structures  such  as  aircraft  skins.  The  probe  contains  a  racetrack 
coil  with  semi-circular  bends  and  linear  straights,  Fig.  1.  The  magnetic  field  between 
the  straights  is  measured  using  a  linear  array  of  magnetic  field  sensors.  The  sensor  array 
samples  the  magnetic  field  at  multiple  sites  without  moving  the  probe  and  hence  reduces 
the  inspection  time.  The  overall  objective  of  this  work  is  to  evaluate  the  capabilities  of 
array  probes  and  assess  their  performance  for  the  detection  of  cracks,  material  loss  and 
surface  roughness  due  to  corrosion.  Here  we  focus  on  the  details  of  the  coil  field  calculation. 

The  theory  for  computing  the  electromagnetic  field  of  a  racetrack  coil,  Fig.  1,  has  been 
developed  by  determining  separately  the  electric  field  due  to  the  bends  and  the  straights 
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It  is  now  possible  to  obtain  ip'  by  substituting  equations  (3.19)  and  (3.21)  into  (3.6). 
Integration  with  respect  to  //,  cf>'  and  z'  gives  a  summation  of  integrals  with  respect  to  k: 


11  °°  4 

V/(r)  =  nl  -  Y,  9X  ,  1  sin[(2A  +  1  )<f>] 

L7r^o2A+1 

roo  ^ 

X  /  [J2\+i(^p)^r2X+i{a,  6,  k)  9  ( k ,  z,  h)  sinh(Kc)]  ck 

Jo 

r° o  ^ 

+  /  [Jo(k/o)^o(«i  b,  k)  9  (k,  z,  h)  sinh(Kc)]  dn 

Jo 


where 


with  ff)  given  by  (3.22)  and 


Jl\z)  =  f  xnJv(x)dx. 

Jo 


These  functions  are  evaluated  for  v  >  2  with  the  aid  of  a  recursion  relationship 


(3.23) 


F„(a,b,K)  =  f  fD{p)Ju{^p)pdp 

Jo 

=  \  ajy\na)  -  bjy\nb)  -  \  J^^ko)  -  jj)2\Kb)  (3.24) 

k  L  J  k ^  L 


(3.25) 


(v  -  n)J™+1(z)  =  —2vznJu(z)  +  (z/  +  n)J„_x{z)  (3.26) 

derived  using  Eq  11.3.6  of  reference  [2]. 

The  integrals  with  respect  to  k  must  be  computed  numerically  and  the  summation  in 
(3.23)  truncated  at  a  suitable  order  depending  on  the  required  accuracy  of  the  result.  For 
the  double-D  filament  loop  [3]  a  truncated  series  of  five  terms  is  sufficiently  accurate  in 
most  cases  and  the  same  is  true  for  the  series  in  (3.23)  representing  the  potential  due  to  a 
racetrack  coil  bend. 

From  equation  (3.2),  the  electric  field  in  cylindrical  coordinates  is 

E(r)  =  -jo W  (p^  -  #£)  ^'(r).  (3.27) 

The  components  of  the  electric  field  are  therefore, 

E  =  _4 gujuonl  ^  cos[(2A  +  1  )(p\ 

nP  t^o 

x  /  J2A+i(K/o)^r2A+i(a,  b,  n)g(n,  z,  zo, /r)  sinh(Kc)  ok  (3.28) 

Jo 

and 


f  1  ^  4 

=  jLup0nI<  -  TTTT  sin[(2A  +  l)4>\ 
l7r  ^o2A  +  1 

r°°  r  2A  + 1  l 

/  kJ2\(kp) - J2\+i(np)  JF2\+i{a,  b,  k)  g  (k,  z,  h)  sinh(Kc)  ok 

Jo  L  P  J 


K,J\{np)JFo{a,  b ,  k)  g  ( k ,  z,  h)  sinh(Kc)  ok 


(3.29) 


respectively. 
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Table  3.2:  Test  Parameters  for  magnetic  sensor  measurements  on  a  hidden  surface  material- 
loss  specimen. 


Coil 


Outer  radius 

10.625  mm 

Inner  radius 

1.6875  mm 

Axial  length 

4.98  mm 

Nominal  lift-off 

2.5825  mm 

Number  of  turns 

337  ±  1  mm 

Number  of  sensors 

33 

Height  of  sensors 

0.869  mm 

Distance  between  sensors 

2.0  mm 

Frequency 

2000  Hz 

Plate 

Conductivity 

1.82  x  107  S/m 

Thickness 

4.85  mm 

Flaw 


Length 

Width 

Depth 


25.4  mm 
25.4  mm 
3.00  mm 


Racetrack  Coil 


31 


FIG.  4.  Magnitude  of  the  magnetic  field  at  32  sensor  sites  due  to  racetrack  coil  excitation 
of  a  metal  plate  containing  a  back  surface  recess,  Fig.  2.  The  excitation  frequency  is  2000 
Hz,  2 d  =  64  mm  (see  Fig.  1)  and  the  other  probe  dimensions  are  as  given  in  Table  I. 
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4.3.2  Dyadic  Form 

Using  (4.1),  (4.8)  and  (4.9)  the  dyadic  Green’s  function  for  the  half-space  conductor  is 
written  as 

£(r|r')  = 

where  r"  =  r'  —  2 zz  is  the  image  point  and 


j+pvv 


G0(r|r')  + 


J/  -  ¥vv' 


G0(r|r")  +  (V  x  z)(V'  x  z)V"( r|r')  (4.16) 


Go(r|r') 


eik  |r-r'| 

47r|r  —  r' 


is  the  familiar  free  space  scalar  Green’s  function  for  the  Helmholtz  equation. 


(4.17) 


4.4  Conducting  Slab 


4.4.1  Scalar  Forms 

For  a  conducting  slab  thickness  d  we  again  adopt  the  form  given  in  equation  (4.11)  for 
both  TE  and  TM  modes  where  we  now  have 


e^(z-z'-2d)  _|_  e-^{z-z  +2d)  I  . 

(4-18) 

with  a  reminder  that  the  chosen  roots  of  y2  =  k2  —  k 2  and  k2  =  u2  +  v2  have  positive  real 
parts.  For  a  conducting  slab  in  air,  reflection  from  top  and  bottom  surfaces  involves  the 
same  reflection  coefficient,  hence,  with  the  upper  sign  for  the  TE  mode  and  the  lower  sign 
for  the  TM  mode 


Gr(z\z')  =  7 -  { e"712"2'1  +  +  g7 {z-z’-id)  + 

2v  l 


Gk{z\z')  =  2-|±e-7(2+2,)  ±e7(2+2'-2d)  +  A(k)  [e7(*+*'-2«0  +  e-7(*+*')J  | 


where 


T(«) 

A(k) 


r2 

1  -  T2e~2~id  ~  1 

r 

i_r2e-27rf 


(4.19) 


(4.20) 


Note  that  once  again 

F'  =  -1  and  T"  =  -  (4.21) 

7  +  K 

for  a  nonmagnetic  slab. 

In  (4.18)  the  second  and  third  term  represent  the  field  due  to  images  created  by  re¬ 
flection  at  the  upper  and  lower  surface  of  the  slab.  Similarly  the  first  two  terms  in  (4.19) 
represent  the  field  due  t  single  reflections  at  upper  and  lower  surfaces.  The  image  terms 
are  separated  in  order  that  the  computation  of  matrix  element  may  be  performed  on  these 
terms  using  the  same  code  that  is  used  for  the  the  contribution  from  the  free  space  Green’s 
function.  This  avoids  taking  them  into  account  via  the  alternative,  Fourier-Bessel  integral, 
representation. 
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2.  Mathematical  model 


Let  us  consider  a  conductive  slab  made  of  a  linear,  isotropic,  non- dispersive,  non¬ 
magnetic,  non- dielectric  and  homogenous  material  of  conductivity  Oo  hosted  in 
free  space.  Let  us  assume  that  the  flaw  is  a  volumetric  flaw  occupying  the  three 
dimensional  domain  Qf  and  that  the  conductivity  of  the  flaw  assumes  the 

constant  value  (see  figure  2.1). 

The  governing  equations  are  (assuming  the  e~jm  time-dependence): 

VxE=  /'(op0H 

V  x  H  =  J  0  +  o  (r  )E  -  j  O)£0E 

V  -E  =  —  (21) 

eo 

V  •  H  =  0 

together  with  the  radiation  condition  at  infinity  and 1 


o 


o 


(r)_  j0°(r)in9*3  vn/ 

V/in^/ 

ja0inflo 

join  9t3  \ 


(2.2) 

(2.3) 


where  is  the  domain  occupied  by  the  conductor  without  the  flaw. 


ECT  coil 


Figure  2.1.  The  reference  geometry.  Q0  is  the  spatial  region  containing  a  homogeneous 
conductive  slab  of  conductivity  G0.  A  homogenous  anomaly  of  conductivity  Of  occupies  the 
volumetric  domain  Cl  f  . 


1  Notice  that  the  symbol  o  0  (r)  stands  for  the  function  defined  by  (2.3)  whereas  o0  stands  for  the 
constant. 
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By  introducing  the  dipole  density  vector  field  P  defined  as: 

P(r)=(a(r)-a0(r))E(r),  (2.4) 

Ampere’s  law  can  be  written  as 

VxH  =  J() +P+o0(r)E- j'co£()E  (2.5) 

showing  that  P  acts  as  a  secondary  source  of  the  electromagnetic  field  produced 
without  the  flaw  (i.e.  £2  =  ). 

Notice  that  P  is  different  from  zero  only  in  £2  (a  (r)  =  o0  (r  )in  91 3  \  £2  f ),  the 
anomalous  region.  As  is  well  known  (see  [1])  P  in  Clf  can  be  found  by  solving 
the  following  integral  equation: 

P(r)  =  P(,)(r)+(a/  -a0)jcop0  jGfe(r  I  r’)P(r')dr',  Vre  Qf  (2.6) 

“/ 

where  P(^  =  (a f  -a0  )e('^  ,  E(,)  is  the  electric  field  due  to  J0  when  a  0  =  a  ,  (i.e. 

without  the  flaw)  and  Gee  is  the  electric-electric  dyadic  Green’s  function  for 
o  f  -  o0 .  Notice  that  G  ee  is  the  solution  of: 

VxVxG“(rl  r')-  jcop0o0 (r)Gee(r  I  r')  =  18 (r - r')  (2.7) 


together  with  the  radiation  condition  at  infinity  and  I  =  xx  +  yy  +  zz  is  the  identity 
operator  expressed  in  terms  of  unit  vectors. 

The  integral  equation  (2.6)  is  interesting  because  its  numerical  solution 
requires  the  discretization  of  £lf  only. 

The  property  underlying  the  numerical  scheme  is  that  P  is  solenoidal  in  Qf  , 
the  interior  of  Q  f .  In  fact,  by  assuming  that  there  is  no  impressed  charge  density 

in  Q  ,  i.e.  V  J()=0  in  £2/  ,  from  Ampere’s  law  (2.1)  it  follows  that 

V  •  E  =  0  in  £2/  (2.8) 


thus  from  (2.5) 

V  •  P  =  (ay  -a0)V-E 
=  0  in  £2  /  . 


(2.9) 


Notice  that  for  given  g0,g  ,  e  (0,+°°)  nothing  can  be  said  on  the  normal  and 
tangential  components  of  P  on  d£2  f . 

The  solenoidality  of  P  in  £2/  can  be  automatically  taken  into  account  by 
introducing  the  dipole  vector  potential  U  defined  by: 
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P  =  VxU. 


(2.10) 


Notice  that  (2.9)  specifies  only  the  curl  of  U  therefore  the  uniqueness  of  U  has  to 
be  imposed  by  gauging  U  as  shown  in  Appendix  B. 


Performance  Report  UTC/AFML  (Project  #  404-25-24)  -  Version  27.9.2002  -  page  44 


where  G„e  is  the  electric-electric  Green’s  function  when  the  whole  space  is  filled 
with  a  conductive  material  of  conductivity  a  0  and  8G ee  is  a  continuous  function 
in  Q0 .  As  is  well  known  [2],G  “  is: 


Goe(rlr')  = 


I  +  -^VV 


g, 


i  (r  I  r') 


(3.8) 


v  J 

where  k 2  =0)2p0(e0  +  jo0/c o )  and  g 0  is  the  scalar  Green’ s  function: 


g 


,(r  I  r)  = 


47t  r-rl 


Following  (3.7)  the  matrix  L  is  the  sum  of  two  terms: 


L  =  Lo  +5L 


where 

-  1*0 1  I v X N » W •  G? (r  I  r) •  V x  N ; (r')drdr' 

af  af 

H  J  J v X  N,  (r). 5G"  (r  I  r1).  V  x N, (r'Jdrdr1 . 

ar  af 


(3.9) 


(3.10) 


(3.11) 

(3.12) 


The  computation  of  8L  can  be  carried  out  by  standard  numerical  integration 
techniques  (the  kernel  SGPC  is  non-singular),  whereas  the  computation  of  L 
requires  volume  and  surface  integrals  (on  dD.  f )  with  kernel  g0  (see  appendix  C): 


L  =  LV+LS 
=0  =0  =0 


(3.13) 


(l =  M  FxNt  (r)*  g0(r  I  r')VxN^(r')drdr'  (3.14) 

-fV  % 

fcol,=-p-  1  Jfi(r).VxN,(r)*0(rlr')nM.VxN;(r')(LS(r')l5(r) 


f  do.  f 


(3.15) 


Notice  that  the  computation  of  the  singular  integral  involved  in  (3.14)  and  (3.15) 
can  be  performed  by  noting  that: 


go 


(r  I  r') 


jk  y*|r-r'| 

—  e 

471 


( 

n  . 
sine 


r 


V 


\ 

+ 


1 

47t  |r  —  r'l 


,  .  /  x  „  sin(x) 

where  sine  (x )  = - 

x 


(3.16) 
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Appendix  A:  Edge  element  shape  functions 


In  the  finite  elements  approach,  the  continuum  is  divided  into  a  finite  number  of 
parts  (elements)  with  a  polyhedral  shape;  every  polyhedral  is  identified  by  the 
coordinates  of  its  vertices  (the  nodes  of  the  discretization).  The  unknown  function 
/( r)  is  represented  in  each  element  by  a  polynomial  approximation  in  terms  of  NP 
parameters  (the  degrees  of  freedom  of  the  solution)  and  of  given  shape  functions: 

/( r)=  £a,W,.(  r).  (A.l) 

i'=l,  AH’ 


Usually,  cij  represents  the  nodal  value  of  /  at  the  z-th  node  of  coordinates  r, 
(a,=/( rf));  each  shape  function  /V,  (nodal  function)  is  locally  based  and  different 
from  zero  only  in  the  elements  sharing  the  same  node  z  and  NP  is  the  number  of 
nodes  in  the  finite  elements  mesh.  In  this  way,  inter-element  continuity  can  be 
assured. 

In  the  following,  we  consider  linear  8-node  isoparametric  brick  elements. 
We  recall  that  the  nodal  function  Nk  associated  to  the  k-th  node  is  continuous, 
piecewise  trilinear,  with  the  following  properties: 

Nk  (r* )  =  1 

Nk  (rm )  =  0,  k  #  m  (A.2) 

I>.(r)  =  l,  tgU, 

j=l,NP 


where  Vj  is  the  finite  element  discretization  of  the  domain  V.  Notice  that  from 
(A.2)  it  follows  that  the  nodal  functions  Nfs  are  linearly  independent. 

We  recall  that  in  this  case  the  /-th  element  is  the  image  of  the  cube  [-1  l]x[-l 
1  ]x[- 1  1]  under  the  mapping  M  (see  Fig.  1)  defined  by: 

kel/ 

(A3) 

keli 

keli 


where  x,y,z  are  the  global  Cartesian  co-ordinates  of  a  point  Q  corresponding  to  the 
point  0  of  local  £,  ,r|,£  co-ordinates  under  the  mapping  M,  I )  is  the  set  of  (eight) 
nodes  of  the  l- th  element,  and  xk,yk,Zk  are  the  global  coordinates  of  the  k- th  node. 
Nk  represents  the  scalar  shape  function  associated  with  the  k- th  node1: 


1  Notice  that  we  use  the  same  symbol  for  the  shape  function  in  the  local  ) ,  local  (rh  sk,  Zp 

and  global  (x,yj)  co-ordinate  systems. 
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system  is  compatible,  the  solution  is  not  unique.  In  fact,  from  (A.  10)  and  (A.  12)  it 
turns  out  that: 

0  =  VxVArt=  £(cg)  S,  (B.7) 

f=lF 

thus1, 


CG  =  0 


(B.8) 


because  the  shape  functions  S  f  are  linearly  independent  (see  Appendix  A).  From 

(B.8)  it  follows  that  adding  any  linear  combination  of  columns  of  G  to  u  is  still  a 

solution  of  (B.6);  in  fact,  the  columns  of  G  yield  irrotational  vector  potential  (and 

zero  dipole  densities)  [3].  Thus,  in  this  case,  imposing  the  gauge  condition  is 
equivalent  to  the  problem  of  assuring  uniqueness  for  system  (B.6). 

Let  us  consider  the  oriented  graph  formed  by  the  nodes  and  oriented  edges 
of  the  edge  element  mesh.  G  and  C  are  the  edge- node  and  face-edge  incidence 

matrices.  From  graph  theory,  the  rank  of  the  ExNP  matrix  G  is  NP- 1  while  the 
rank  of  the  FxE  matrix  C  is  E-NP+ 1.  The  uniqueness  of  the  solution  of  (B.6)  is 

guaranteed  if  U  in  (B.2)  is  represented  by  E-NP+ 1  degrees  of  freedom. 

To  select  the  proper  degrees  of  freedom,  we  notice  that  the  uniqueness  of  the 
solution  of  (B.6)  means  that  VxU^  0  when  at  least  one  of  the  E-NP+l  degrees 
of  freedom  is  different  from  zero.  Let  us  decompose  the  graph  into  an  arbitrary 
tree  (NP- 1  edges)  and  cotree  (the  residual  E-NP+l  edges)  (the  algorithm  to 
identify  all  possible  trees  of  a  graph  is  presented  in  [6]).  Then,  assuming  Ue= 0 
when  e  is  an  edge  of  the  tree,  we  notice  that  V  x  U  ^  0  if  at  least  one  degree  of 
freedom  related  to  a  cotree  edge  is  different  from  zero.  In  fact,  each  edge  of  the 
cotree  closes  a  single  independent  loop  y  e  with  the  edges  of  the  tree  and  the 
circulation  of  U  along  ye  is  equal  to  ±Ue  (y ,,  does  not  include  cotree  edges  apart 

from  edge  e).  Therefore,  the  E-NP+l  degrees  of  freedom  associated  to  the  edges 
of  the  cotree  lead  to  a  reduced  set  of  equations  that,  if  compatible,  admit  a  unique 
solution  [5,  7,  8]. 

For  the  sake  of  completeness,  it  is  worth  noting  that  the  NP-l  degrees  of 
freedom  associated  to  the  tree  edges  can  be  used  to  represent  gradient  vector 
fields.  In  fact,  once  the  degrees  of  freedom  of  the  edges  of  the  tree  have  been 
assigned,  for  any  cotree  edge  e  we  choose  Ue  so  that  the  circulation  of  U  along  y  e 
is  equal  to  zero  (since  the  circulation  of  U  along  y e  involves  only  the  cotree  edge 
e ,  we  find  that  Ue  depends  linearly  on  the  degrees  of  freedom  associated  to  the 
tree  edges).  Therefore,  (B.2)  gives  a  vector  field  U  that  is  irrotational. 


1  Notice  that  (B.8)  is  the  discrete  equivalent  of  VxV/=0. 
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Assembling  V,  E,  L,  8L,  L  and  L  SQ 

The  integrals  appearing  in  the  numerical  formulation  can  be  written  in  the 
following  way: 


(v'T =£>„, 

m= 1 

(1) 

Ne 

(E)»;=Ep.,.J 

m= 1 

(2) 

N  N 

^)kj  ~ 

m= 1  l—\ 

(3) 

where  Ne  is  the  number  of  elements  of  the  mesh, 
vm,,  =  JVxN,(r)-E(0(r)dr 

T,„ 

(4) 

P m.k,i  = — 1 — JvxN((r)  VxN;(r)dr 

rr  — rr  J 


X 


(5) 


a 


p0  J  J  Vx  N,  (r  )•  Gee{r  I  r’)  •  Vx  N  .  (r')dT  dr  . 


(6) 


and  x,„  is  the  three  dimensional  domain  of  the  m-th  element  of  the  finite  element 
mesh. 


The  standard  procedure  for  assembling  V,  E  and  L,  assuming  that  the  non¬ 
singular  integrals  (4)  and  (5)  are  computed  by  the  Gauss  quadrature  method,  is 
given  by  the  following  algorithms  ( NG  is  the  number  of  Gauus  points  per  element, 

I^E  is  the  set  of  indices  of  the  cotree  edges  lying  on  the  boundary  of  element  m ): 


Procedure  1  :  compute  and  store  the  values  of  the  shape  functions  assumed  in  the 
Gauss  points  for  any  element 


for  m—1,  Ne 


for  ke  I 


CE 

m 


for  n=l,  ...,Ng 

set  SHPF(m,k,n)=  V x N, (r„m ) 


1  For  given  m,  k  and  n,  SHPF(m,k,n)  is  a  three  component  vector. 
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SamM,,=|i0jjVxN^r)-5Ge1rlrO-VxN.(rVdr, 

tmX, 

(7) 

whereas  in  the  second  case  amJ  k  .  has  to  be  replaced  by  CL^l  k  j  defined  as 
=  |i.JJVxNiW-*o(r  I  r')VxN;(r')dr'dr 

xmx, 

(8) 

The  numerical  computation  of  8am ,  k  .  can  be  carried  out  by  Gauss  quadrature 
rule,  whereas  the  numerical  computation  of  CL°lkj  can  be  carried  out  by  taking 

into  account  (3.16).  Specifically,  the  term  arising  from  r^“sinc 

4jc 

can  be  computed  by  the  Gauss  quadrature  rule,  whereas  the  term  arising  from 

— - T  can  be  computed  by  the  approach  proposed  in  [8]. 

4n  |r  -  r'| 

Finally,  the  surface  integral  (3.15)  can  be  written  as: 

/  \  %  Nhf 

i'-:l  -XIX, 

m=  1  /  =1 

(9) 

where  A^/  is  the  number  of  boundary  faces  of  the  mesh, 

P,,r.,j--7rjl"(r) 'VxN,(r)s,(r  'r')ft(r') .VxN;(r')dS(r')lS(r) 

* 

(10) 

and  Em  is  the  m-th  boundary  face.  The  algorithm  for  computing  is  (I~E  is  the 
set  of  indices  of  cotree  edges  lying  on  the  boundary  surface  in  ): 

Procedure  5  :  compute  L*  1 
set  L  |  =0; 

compute  and  store  the  values  of  n(r)- VxNA  (r)  evaluated  at  the  Gauss 
points  of  any  boundary  face; 
for  m  =1,  . . .,  Nbf 

1  The  algorithm  for  computing  p  _  ,  .  is  left  to  the  reader  (for  instance,  combine  (3. 16)  with  (18) 
and  (19)  in  [8]). 
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The  incidence  matrices  and  other  relevant  geometrical  and  topological  data 


Input 

x(i,  1:3),  z=l  :NP 
ix(i, i=l:S,j=l:Ne 

ix(9,  /'),  7=  1 : 

Output 

ixb(z,j),  i  =VA,j=l:Nbf 
ltact(z,y),  z'=l:12,j=l:./Ve 

ltact_bou(z,  /),  z=  1 :4, y=  1 

inc  ( 1 ,  j) ,  inc(2 ,  j) ,  j=  1 

ncoalb(z),  z=l:iVC£ 
latel(z,j),  i  =l:l2,j=l:Ne 

ielface  (j) ,  j=  1  :Nbf 

indfac e(z,j),  i=l:4,  j=l:  A/*/ 


coordinates  of  the  z-th  nodal  point 
label  of  the  z-th  node  of  the  j-th 
element  (see  fig  2) 

label  of  the  material  properties  of  the 
j-th  element 


label  of  the  z-th  node  of  the  / - th 
boundary  face 

label  of  the  z-th  cotree  edge  of  the  j- 

th  element,  ltact(zj)  is  zero  if  the  z-th 

edge  belongs  to  the  tree. 

label  of  the  z-th  cotree  edge  of  the  j- 

th  boundary  face.  ltact_bou(z,  j)  is 

zero  if  the  z-th  edge  belongs  to  the 

tree. 

labels  of  the  two  nodes  identifying 
the  j-th  edge 

label  of  the  z-th  edge  of  the  co-tree 
label  of  the  z-th  edge  of  the  j-th 
element  (see  fig  2) 

label  of  the  element  of  the  j-th 
boundary  face 

(local)  label  of  the  z'-th  edge  of  the  j- 
th  boundary  face  (see  fig.  2) 


Procedure  1  :  find  boundary  faces  (Nbj,  ixb,  ielface  and  indface) 
j=0; 

for  iell=l,  ...,Ne 

for  each  face  jp  of  element  iM 
set  ind=false; 
while  ind=false 

for  iej2=  ieii+1,  ■■■,  Ne 

for  each  face  jp  of  element  iei2 
if  jfi=jj2  set  ind=true;1 

if  ind=false 

j=j+i; 

insert  the  nodes  of  face  jp  in  the  listixb(l:4,  j); 
insert  i eii  in  the  listielface(j); 
insert  jp  in  the  list  indface  (1:4,  j); 

Nbj=ji 


'  jfi  =jfi  means  that  the  set  of  nodes  of  face  jf,  is  equal  to  the  set  of  nodes  of  face  jp. 
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Summary 

The  central  task  of  Phase  II  is  to  reconfigure  the  numerical  scheme 
using  edge  elements.  The  motivation  for  this  is  that  the  array  probe 
simulation  is  then  based  on  an  approach  that  is  more  reliable  and  efficient 
than  the  previous  method  based  on  traditional  volume  elements.  With 
an  improved  algorithm,  calculation  of  the  electromagnetic  field  in  such 
structure  as  a  fastener  or  a  lap  joint  can  be  performed  faster  because 
the  underlying  computational  engine  uses  as  few  unknown  as  possible 
and  is  more  reliable  because  it  is  constructed  in  such  a  way  as  to  avoid 
spurious  modes. 

In  fact  the  numerical  approach  that  represents  the  field  in  terms  of  edge 
elements,  was  develop  over  the  last  decade  to  overcome  the  problem 
of  spurious  solutions.  The  group  at  Cassino,  who  participated  in  the 
present  project,  were  at  the  forefront  of  this  effort  and  we  are  fortunate 
that  they  could  partner  with  CNDE  to  create  a  computational  scheme 
that,  when  fully  developed,  will  be  at  the  forefront  of  modern  numerical 
techniques. 

Due  to  delays  in  finalizing  the  agreement  between  Iowa  State  University 
and  Cassino,  the  work  on  the  new  numerical  formulation  was  not  com¬ 
pleted  until  August  which  left  little  time  for  coding.  However,  Vipul 
Katyal  has  work  hard  on  the  prototype  of  the  new  code  at  ISU,  and 
initial  tests  for  bugs  accuracy  and  performance  have  been  carried  out. 
Further  work  is  needed  before  the  code  is  validated  and  can  be  used  for 
routine  calculations. 

Improvement  in  the  experimental  methods  for  validation  have  been  made 
by  using  a  Stanford  lock-in  analyzer  for  magnetic  sensors  measurements. 
New  probes  have  been  built  for  the  Phase  II  validation  work  with  an 
improve  signal  to  noise  performance.  Because  the  measuring  systems  has 
been  upgraded,  new  measurements  have  been  performed  on  specimens 
manufactured  in  the  first  year  of  this  project  and  some  new  test  pieces 
added. 

This  report  is  divided  into  a  series  of  articles  beginning  with  a  summary 
of  some  numerical  and  experimental  results.  This  is  followed  by  a  brief 
account  of  a  proposed  imaging  method  for  assessing  material  loss  and 
Article  3  is  a  copy  of  a  paper  on  the  racetrack  coil  submitted  to  Jour¬ 
nal  of  Applied  Physics.  In  addition  we  have  included  a  technical  review 
of  dyadic  Green’s  functions  for  layered  structures,  Article  4,  which  was 
written  to  coordinate  theoretical  developments  at  Cassino  and  ISU.  Fi¬ 
nally  we  include  a  comprehensive  report  from  Rubinacci  and  Taburrino 
giving  the  edge  element  theory  and  an  outline  of  the  new  numerical 
formulation. 


VI 
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Figure  3:  Asahi  Kasai  magnetic  field  sensor.  Dimensions  are  given  in  millimeters. 


sensor  current,  the  result  is  largely  independent  of  inductive  pick-up  signals. 

1.2.1  Calibration  Test 

Two  precision  probes,  whose  parameters  are  given  in  Table  1,  have  been  made  and  charac¬ 
terized  following  a  procedure  given  in  Appendix  A.  The  characterization  ensures  that  the 
probe  performance  is  in  accordance  with  the  assumptions  of  the  simulation  and  that  the 
sensors  is  calibrated. 

In  the  characterization  procedure,  a  test  is  carried  out  to  ensures  that  the  value  adopted 
for  the  coil  probe  liftoff  yields  predictions  of  coil  impedance  and  its  variation  with  frequency 
that  are  consistent  with  measurements  made  with  the  probe  on  a  homogeneous  copper  slab 
(See  Appendix  A  for  these  results).  An  additional  test  is  performed  to  ensures  that  the 
value  adopted  for  the  Hall  probe  liftoff  yields  theoretical  predictions  of  the  magnetic  field 
variation  with  frequency  at  the  sensor  site  that  are  consistent  with  experimental  field 
measurements  for  a  probe  on  the  same  unflawed  copper  slab. 

In  this  latter  test,  the  axial  magnetic  field  at  the  sensor  has  been  predicted  for  a 
probe  on  a  copper  half-space  and  compared  with  calibrated  measurements  made  with  the 
probe  an  a  thick  copper  plate.  One  assumes  that,  because  the  theory  is  fairly  basic, 
the  theoretical  predictions  are  correct.  Then  the  experimental  data  is  examined  to  see  if 
the  calibration  is  accurate  and  the  measurement  system  is  functioning  correctly  without 
introducing  spurious  effects.  The  results  of  this  comparison,  Figure  4,  show  that  the 
experimental  field  measurements  and  predictions  agree  to  with  4%  over  the  operating 
frequency  range  of  the  probe. 

1.2.2  Square  recess 

The  basic  test  case  for  the  determination  of  material  loss  is  a  square  recess  in  the  bottom 
surface  of  a  aluminum  plate.  A  comparison  between  theory  and  experiment  for  this  case 
was  carried  out  in  Phase  I  and  gave  reasonable  results.  However  there  was  a  tendency  for 
the  model  to  overestimate  the  magnitude  experimental  field  measurements.  With  the  aim 
of  seeking  improved  agreement  between  theory  and  experiment,  new  data  was  obtained  on 
the  square  recess  specimen  using  a  new  probe  and  an  improved  calibration  procedure.  In 
addition,  the  analysis  of  the  field  predictions  was  reviewed  and  modified. 

The  recess  is  25.45  mm  square  in  the  lower  face  a  test  piece  4.85  mm  thick.  The 
recess  was  3  mm  deep  which  means  that  the  material  loss  is  detected  through  1.85  mm 
of  aluminum.  Phase  and  amplitude  measurements  of  the  magnetic  field  at  the  site  of  the 
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Figure  5:  Magnetic  field  measurements  of  the  signals  due  to  a  square  recess  at  a  frequency 
of  2.5kHz  using  probe  P3. 


plate,  4.8  mm  thick,  is  shown  in  Figure  7.  These  measurements  will  be  compared  with 
theoretical  predictions  in  the  future.  They  were  made  in  order  create  a  stepping  stone 
towards  the  problem  of  computing  and  validating  the  response  for  a  fastener  in  a  plate, 
Figure  8.  However,  time  permitted  an  examination  of  the  field  predictions  for  a  fastener 
only,  Figure  8. 

The  calculations  were  performed  by  defining  a  region  whose  conductivity  differed  from 
that  of  the  host  aluminum  plate.  Because  it  proved  difficult  to  measure  the  conductivity  of 
the  fastener,  results  were  calculated  using  a  number  of  different  fastener  conductivities  and 
the  best  results  chosen,  Figure  9.  This  provides  a  test  of  the  predicted  shape  rather  than 
magnitude  since  the  latter  was  adjust  through  changes  of  fastener  conductivity  to  give  the 
best  fit. 

1.2.4  Lap  joint 

A  basic  lap  joint  structure,  Figure  10  was  used  to  simulate  the  P3  probe  signal  at  500Hz 
and  2500  Hz  as  shown  in  Figure  11.  These  results  show  firstly  that  the  signal  has  a 
strong  frequency  dependance  and  secondly  that  the  higher  frequency  signal  gives  a  sharper 
edge  response  and  therefore  may  be  said  to  provide  a  better  resolution  of  the  underlying 
structure. 

A  comparison  between  theoretical  predictions  of  the  simulation  and  measurements  is 
shown  in  Figure  12.  The  magnitude  of  the  predictions  are  within  about  15%  of  the  mea¬ 
surements.  There  is  a  small  but  significant  phase  error  and  some  detailed  features  of  the 
measurements  do  not  appear  in  the  theoretical  predictions.  This  may  be  due  to  the  fact, 
common  to  almost  all  numerical  models  at  this  time,  that  sharp  edges  give  rise  to  edge 
singularities  and  these  are  not  properly  accounted  for  in  discrete  representation  of  the  field 
used  for  computation. 
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Figure  8:  Model  representation  of  a  5.5  mm  diameter  countersunk  rivet  in  a  conducting 
plate. 


Figure  9:  Magnetic  field  measurements  of  a  5.5  mm  diameter  rivet  in  a  conducting  plate, 
comparison  of  experiment  with  theory.  The  unknown  fastener  conductivity  was  varied  for 
the  best  fit. 


Figure  10:  Simplified  lap  joint:  2mm  thick  plates,  45  mm  wide  lap  joint  with  a  0.14  mm 
insulating  layer  between. 
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Appendix  A 

Probe  Characterization  and  Calibration  Procedure 

1.  Coil  winding:  The  coil  former  dimensions  are  measured,  the  coil  is  wound  evenly 
on  the  former,  the  number  of  turns  noted  and  the  outer  diameter  measured. 

2.  Free  space  impedance:  The  coil  impedance  is  measured  in  free  space  as  a  function 
of  frequency  over  the  operating  frequency  range  of  the  probe  (40  Hz-10,000  Hz  in  this 
case)  using  an  impedance  analyzer  (Agilent  4294A).  The  free  space  self  inductance  of 
the  probe,  Lq  in  the  zero  frequency  limit  is  found  by  extrapolating  the  self-inductance 
data  to  the  zero  frequency  limit.  The  theoretic  value  of  Lq  is  compared  with  the 
experimental  value  to  check  that  it  is  consistent  with  the  coil  dimensions  and  the 
number  of  turns.  If  the  calculated  Lq  disagrees  with  measurement,  the  outer  radius 
of  the  coil  is  incrementally  adjusted,  Figure  A1  until  agreement  is  reached. 

3.  Correction  inter  winding  capacitance  The  effective  shunt  capacitance  of  the  coil, 
Co  is  found.  This  value  is  used  to  correct  measured  impedance  before  comparison 
with  theory. 

4.  Half-space  coil  impedance:  The  coil  impedance  is  measured  as  a  function  of 
frequency  with  the  probe  on  a  thick  conducting  copper  plate.  Measurements  are 
made  over  the  operating  frequency  range  of  the  probe.  These  measurements  are  then 
corrected  for  the  effect  of  the  shunt  capacitance  and  compared  with  the  theoretical 
predictions  from  the  theory  of  Dodd  and  Deeds  [3].  A  root  mean  square  error  function 
is  used  to  quantifying  the  fit  of  all  the  data  to  the  theory  and  is  calculated  from 

■£2=  £  \ztp,)-zM\2  (Ai) 

m=l,M 

where  ZmXptr>  is  the  corrected  experimental  measurement  at  frequency  u ;m,  Z(com)  is 
the  theoretical  value  and  the  sum  is  over  for  all  M  measurements.  The  effective  coil 
lift-off  value  zc  is  then  adjusted  until  this  error  is  a  minimum. 

5.  Sensor  location:  The  Hall  sensor  is  then  mounted  on  the  coil  axis  and  the  location 
of  the  sensitive  region  found  from  a  micro-focus  X-ray. 

6.  Magnetic  sensor  calibration:  The  Hall  sensor  is  calibrated  in  air  as  follows.  A 
standard  resistor  is  connected  in  series  with  the  coil  and  the  ac  signal  output  from  a 
Stanford  lock-in  analyzer  applied  across  the  coil  and  resistor.  The  lock-in  is  used  to 
determine  the  coil  current  and  the  Hall  voltage.  Because  the  coil  is  in  air,  a  simple 
expression  can  be  used  to  determine  the  magnetic  field  at  the  site  of  the  Hall  sensor. 
Since  the  field  is  known,  the  sensor  can  be  calibrated  over  the  operating  frequency 
range  of  the  probe. 

7.  Half-space  field  measurements:  The  field  measured  by  the  Hall  device  is  then 
determine  as  a  function  of  frequency  with  the  probe  on  a  thick  metal  copper  plate 
of  known  conductivity.  These  results  are  then  compared  with  the  theoretical  values 
and  a  least  squares  error  function  computed  as  for  the  impedance  data: 
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Figure  Bl:  Coil  resistance  and  self-inductance  variation  with  frequency. 

B2  Sensor  Calibration 

The  coil  current  excitation  is  provided  by  a  lock-in  amplifier  output  and  measured  using 
the  lock-in  by  monitoring  the  voltage  drop  across  a  10  H  standard  resistor  in  series  with  the 
coil.  Similarly  the  Hall  sensor  voltage  is  measured  by  the  lock  in.  Because  the  instrument 
acts  as  a  voltage  source  and  the  load  is  inductive,  the  load  impedance  and  the  coil  current 
varies  with  frequency  and  hence  the  Hall  signal  varies  with  frequency  as  shown  by  probe 
measurements  in  free  space,  Figure  B3. 

The  ratio  of  the  Hall  voltage  to  the  coil  current  should  be  reasonably  constant  provided 
that  the  operating  frequency  range  is  within  the  bandwidth  of  the  Hall  sensor  preamplifier. 
This  is  certainly  true  for  the  dominant  real  part1,  Figure  B4.  It  is  estimated  that  the 
Hall  sensor,  like  the  coil,  has  a  lift-off  of  0.71  mm.  Using  this  value  and  the  other  probe 
parameters,  the  field  on  the  coil  axis  at  the  site  of  the  sensor  has  been  calculated.  The 
calculated  value  of  the  field  at  the  sensor  for  the  probe  in  free  space  per  unit  coil  current  is 
Ho  =  81, 316Amps/m.  Note  also  the  characteristic  parameter  Rhc  =  5. 31714  which  is  the 
real  part  of  the  ratio  of  the  Hall  voltage  to  the  probe  current  with  the  probe  in  air.  This 
value  varies  slightly,  see  Figure  B4,  but  we  adopt  a  constant  value  determined  at  2000  Hz. 
From  measurements  of  the  Hall  sensor  voltage  Vh  and  the  coil  current  Ic,  the  magnetic 
field  at  the  sensor  is  given  by 


(B1) 

Using  this  calibration  procedure,  the  axial  magnetic  field  at  the  sensor  has  been  pre¬ 
dicted  for  a  probe  on  a  copper  half-space  and  compared  with  calibrated  measurements 
made  with  the  probe  on  a  thick  copper  plate.  These  results  are  show  in  the  main  text. 


1  The  real  part  corresponds  to  the  component  in  phase  with  the  reference  signal  of  the  lock-in. 
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Figure  B2:  Coil  impedance  variation  with  frequency;  comparison  of  experiment,  circles, 
with  theory,  line. 
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where  the  Green’s  function,  like  ip' ,  is  continuous  at  the  air  conductor  interface,  has  a 
continuous  normal  gradient  and  vanishes  at  infinity.  The  Fourier  transform  with  respect 
to  x  and  y  is  written 

rsj  rco  rco 

f(u,v)=  /  f(x,  y)e~iux~ivy  dx  dy.  (3.7) 

J  —  oo  J —co 

Hence  by  taking  the  Fourier  transform  of  equation  (3.6)  with  respect  to  x  and  y,  and  noting 
the  convolutional  properties  of  the  integral,  it  is  found  that 

ip'  (u,v,  z)  = —  /  9  (u,v,  z,  z')  m(u,v,  z1)  dz' ,  (3.8) 

/io  Jh—c 

where  the  integration  is  between  the  lower  and  upper  limits  of  the  source  coil  and  9  and  m 
are  the  Fourier  transforms  of  G  and  M  respectively. 

The  y-conrponent  of  the  current  density  in  the  straight  elements  of  the  source  coil,  Fig. 
1,  is  written  as 


Mr)  = 


nl sign(x),  h  —  c  <  z  <  h  +  c,  b  <  |x|  <  a,  \y\  <  d, 

0,  otherwise 


(3.9) 


where  I  is  the  current,  n  the  number  of  turns  per  unit  area  and  the  current  is  deemed  to 
flow  in  a  counter-clockwise  direction  viewed  from  above.  It  can  be  deduced  from  (3.3),  by 
writing 

Ho  nlf(x,y),  li  —  c<z<h  +  c 

0,  otherwise 


M(x,y,z)  = 


(3.10) 


that 


f{x,y)  = 


0  <  |x|  <  b,  \y\  <  d, 

b  <  \x\  <  a,  \y\  <  d, 

otherwise. 


(3.11) 


Because  f(x,  y)  is  even  in  x  and  y,  the  Fourier  transform  may  be  written  in  the  form  of 
the  double  cosine  integral 


r  oo  rco 

f(u,v)  =  4  /  /  f  (x ,  y)  cos(ux)  cos(vy)  dx  dy 

Jo  Jo 

4 

= - y~  [cos(ua)  —  cos(u6)]  sin(vd). 


u*v 


(3.12) 


An  expression  for  the  electric  field  which  can  be  evaluated  numerically  is  obtained  in 
the  following  way.  The  Fourier  transform  of  (3.2)  with  respect  to  x  and  y  is 


e  (u,  v,  z)  =  —buy, o(vx  —  uy)  ij)'  ( u ,  v,  z) 
where  is  given  by  (3.8)  and  (3.10)  as 

~  ~  rh-\-c  ^ 

ip1  (u,v,  z)  =  nl  f  (u,v)  /  9  (u,v,  z,  z')  dz' ,  z  <  0. 

J  h—c 


(3.13) 


(3.14) 
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For  a  half-space  conductor 

(3-15) 

7  +  K 

where  7  =  \J k2  +  jug^o,  taking  the  root  with  a  positive  real  part.  Similarly,  Green’s 
function  for  a  slab  is  computed  by  taking  into  account  the  reflection  from  the  internal 
surfaces, 

~  ,  1  , 1  +  Te~2l^d+Z^ 

Ww')  =  l _  r2e-27,  (3-16) 

where  F  =  ■,  the  reflection  term  and  d  =  height  of  the  slab.  Performing  the  integration 

in  (3.14)  gives 

~  ~ 

ip  (u,v,z)  =  2nl  f  (u,v)  g  (u,v,z,h)sinh.(Kc)/K,  z  <  0,  (3-17) 

hence  by  substituting  into  (3.13)  it  is  found  that 

e  (it,  v,  z)  =  —  2ujgonI - f  (u,v)  9  (u,v,  z,h)smh.(nc),  z<  0.  (3.18) 

The  electric  field  can  now  be  computed  using  a  fast-Fourier-transform  algorithm. 


3.3  D-Coil  Field 


Two  D-coils  are  used  to  represent  the  bends  of  the  racetrack  coil,  Fig.  1.  Using  essentially 
the  same  formulation  that  was  used  for  the  linear  coil,  ip'  for  the  D-coil  is  written  as  in 
(3.6).  It  is  convenient  to  express  the  Green’s  function  in  cylindrical  polar  coordinates,  as 


1  OO 

m=  0 


,cos  [m(<p  —  <f>')\  J  Jm{np)Jm{Kp')  g  (re,  z,  z')ndK,  z<  0,  (3.19) 


which  can  be  derived  using  an  approach  given  by  Morse  and  Feshbach  [1],  In  (3.19),  em  is 
the  Neumann  factor:  eo  =  1  and  em  =  2  (m  =  1, 2,  3, . . .).  For  the  interior  of  a  half-space 
conductor,  g  is  given  by  (3.15).  For  for  a  slab,  g  is  given  by  (3.16).  In  order  to  evaluate 
(3.6),  the  explicit  form  of  M(r)  appropriate  for  the  D-coil  must  be  found.  This  form  is 
developed  as  follows. 

The  azimuthal  counter-clockwise  current  in  a  D-coil  is  written  as 


Mr) 


nl ,  h  —  c  <  z  <  h  +  c,  0<</><7T,  b  <  p  <  a 

0,  otherwise. 


(3.20) 


where  /  is  the  current,  n  the  number  of  turns  per  unit  area.  It  can  be  deduced  from  (3.20), 
by  writing 


that 


M(p,  <p,  z) 


Po  nlfnip),  h  —  c<z<h  +  c,  0  <  (p  <  ir, 

0,  otherwise 


(3.21) 


Id(p )  = 


a  —  b, 
a  —  p. 
0, 


0  <  p  <  6, 

b  <  p  <  a, 
otherwise 


(3.22) 


Racetrack  Coil 


FIG.  1.  Racetrack  probe  showing  coil  geometry  and  magnetic  field  sensor  array. 


FIG.  2.  Conducting  plate  with  a  square  recess 
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'  r)C' 

^  =0  and  \nG"]  =  0  (4.4) 

Similarly  the  continuity  of  Hj  implies  that 

T  r)G "  1 

[<jG'}  =  0  and  =  0.  (4.5) 

The  continuity  conditions  are  applied  to  determine  the  Fourier  transform  of  the  TE  and 
TM  Green’s  functions. 

4.2.3  Fourier  Representation 

Using  a  two-dimensional  Fourier  representation,  denoted  by  a  tilde,  we  define 

G(z,z')=  r  [°°  Girlr'y-^-^-^y-^dudv,  (4.6) 

J — OO  J — oo 

with  the  inverse 

1  roc  roc  _ 

G(r\r')  =  -—,  /  G{z,z'Yu{-x-x">+i<y-y,)dudv,  (4.7) 

(2tt J  J — oo  J — oo 

where  the  tilde  implies  dependence  on  the  Fourier  space  coordinates  u  and  v. 

4.2.4  Further  Useful  Dyadic  Forms 

The  scalar  decomposition,  (4.1)  is  transformed  to  more  convenient  forms  by  using  the 
following  two  equations  which  can  be  derived  using  identities  give  on  page  18  of  the  text 
by  Felsen  and  Marcuvitz  [2]. 

1+l 2VV  GT(r|r')  =  ^(Vx\/xz){VxVxz)Ur(r\r')  +  (\/xz){Vxz)Ur{r\r')  (4.8) 

where  the  subscript  T  denotes  ,  x  —  x1,  y  —  y' ,  z  —  z'  dependence.  Although  the  scalar 
form  (4.1)  and  the  equivalence  (4.8)  are  not  valid  at  the  singular  point,  our  final  result 
will  contain  the  correct  representation  of  the  singularity  in  the  form  of  the  left  hand  side 
of  (4.8).  The  second  equation  for  transforming  (4.1)  is 

1'  -  pVV'  GA(r|r")  =  (V  x  V  x  z)(V'  x  V'  x  z)UA{r\r")  +  (V  x  z)(V  x  z)UA(r\r") 

(4.9) 

where  the  subscript  A  denotes  a  x  —  xf,  y  —  y’ ,  z  +  z!  dependence  and  1'  =  xx  +  yy  —  zz. 

4.2.5  Evaluation  of  the  Scalar  Green’s  Function 

Suppose  we  have  a  multilayered  structure  then  we  need  a  subscript  to  identify  which  layer 
is  being  referred  to.  We  shall  not  however  introduce  this  layer  subscript  since  the  details 
will  be  omitted.  However,  to  summarize  briefly,  we  locate  the  singular  source  in  one  of 
these  layers.  Then  proceed  as  follows: 

•  Write  the  Fourier  transform  of  the  solution  of  (4.3)  as 

G{z,z')  =  A^e-^  +  B{K)e^z 


for  each  layer. 
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•  Apply  the  boundary  conditions  to  find  all  the  A  and  B  coefficients.  If  the  number  of 
layers  above  and  below  the  source  layer  is  arbitrary  then  the  relationships  between 
the  coefficients  can  be  expressed  recursively  [5]  as  is  done  for  example,  by  Dodd, 
Cheng  and  Deeds  [6]  for  the  field  in  a  multilayered  cylindrical  structure. 

•  Use  (4.1)  (4.8)  and  (4.9)  to  express  the  dyadic  Green’s  function  in  a  form  suitable 
for  computation. 

Two  examples  are  given  below. 

4.3  Homogeneous  Half-Space 

4.3.1  Scalar  Forms 

A  simple  example  of  the  general  layered  domain  is  one  where  the  source  is  located  in  a 
homogeneous  half-space  conductor  (z’jO)  adjoining  a  non-conducting  half-space  (z^O).  For 
the  field  in  the  conducting  half-space  (z  <  0), 


G'(z,z') 

1 

e-7|z-2'|  _  e7 (z+z’) 

27 

- 

G"(z,  z') 

1 

\-l\z-z' I  ,  7-  K^(z+z')] 

27 

0  0 

7  fK 

J 

with  y2  =  k2  —  k2  and  n2  =  u2  +  v2,  taking  positive  real  roots  in  each  case. 

As  in  more  complicated  cases,  such  as  the  slab  below,  it  is  useful  to  express  the  scalar 
Green’s  functions,  for  both  TE  and  TM  modes,  as  (suppressing  the  primes  distinguishing 
TE  and  TM  modes) 


G{z\z')  =  Gr{z\z')  +  G\{z\z!),  (4.11) 

where,  the  z  —  z'  dependent  Green’s  function  Gy(z\z')  is  just  the  unbounded  domain  term, 
Go  say, 

Gr{z\z')  =  'I  =  G0(z\z'),  (4.12) 

whereas  the  z  +  z1  dependent  term  is 

Ga  =  ^T(K)e^z+*')  =  -G0(z\  -  z')  +  V(z\z')  (4.13) 

27 

where 

V(z\z')  =  ?  !'(/.•)  +  l}e^z+z’)  (4.14) 

27 

In  this  way,  the  modes  are  distinguished  by  having  different  ‘reflection’  coefficients: 

T'  =  -1  and  T"  =  -  (4.15) 

7  +  K 


for  a  nonmagnetic  conductor.  Note  that  because  the  TM  reflection  coefficient  is  —1  V'  =  0. 
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1  p  OO  1  , 

w+(p,z,z')  =  -—  — A'(k)  e7(*+*-2«0  +  e-7(*+*)  j0(K/0)K£k  (4.32) 

Z7T  J 0  Z7  L  J 

1  /*oo  1  r  /  i 

F_(p,  z,z')  =  —  - [T"(/e)  -  T'(/e)l  \e^z~z  +  e-rtg-z  +2d)  JQ(np)dK  (4.33) 

27t  Jo  2/^7  L  J  L  J 

1  p  oo  1 

^+(p,  z,z')  =  —  —  [A"(«)  +  A'(k)]  -“)  +  e-^+y)  JqM^,  (4.34) 

Z7T  Jo  ^7  J 

where  p  =  [(x  —  x')2  +  (y  —  y')2]1/2- 


Summary 


1.  Introduction 

2.  Mathematical  model 

3.  Numerical  model 

Appendix  A:  edge  element  shape  functions 
Appendix  B:  Tree-Cotree  Decomposition 
Appendix  C:  integration  of  the  hyper-singular  term 
References 


Performance  Report  UTC/AFML  (Project  #  404-25-24)  -  Version  27.9.2002  -  page  40 


3.  Numerical  model 


The  numerical  model  proposed  in  this  section  is  obtained  by  applying  Galerkin’s 
method  to  eq.  (2.6).  The  unknown  is  the  dipole  vector  potential  U  which  we 
represent  as  the  linear  combination  of  edge  element  based  shape  functions,  i.e. 

U  =  ^t/eNe  where  the  Ne‘s  are  edge  element  based  shape  functions  and  e 

e 

represents  an  edge  of  the  mesh  (see  Appendix  A).  Specifically,  the  sum  is 
restricted  to  the  edges  of  a  cotree  of  the  finite  element  to  numerically  impose  the 
uniqueness  of  U  (see  Appendix  B). 

The  unknown  dipole  density  (see  2.10)  is  represented  as: 

P  =  E'7.P.  (3-D 


where  Pe  is  expressed  in  terms  of  Ne  as: 

Pe=VxN„.  (3.2) 

Galerkin’s  method  applied  to  (2.6)  gives: 

(R-  jcoL)u  =  V(,)  (3.3) 

where  U  is  the  column  vector  made  of  the  Ue‘ s  and1 

(v'T  =  |vxN,.(r).E(l,(r)dr 

A 

VxNj(r)'VxN  j  (r  )dr 
J  JVxN1(r).G"(rlr')  VxNJ(r')dr'dr 

af  A 

Notice  that  the  computation  of  the  matrix  L  can  be  troublesome  due  to  the 
hyper- singular  part  of  the  dyadic  Green’s  function  Gee.  However,  the  computation 
of  hyper-singular  integrals  can  be  reduced  to  the  computation  of  surface  integrals 
on  dQ  f  .  In  order  to  show  this,  we  first  notice  that 

GeP(r  lr’)=G"(r  I  r')  +  5Gee(r  I  r')  (3.7) 


(3.4) 

(3.5) 

(3.6) 


1  In  (3.4)-(3.6)  k  and  j  represent  edges  of  the  finite  elements  mesh. 
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(A.  4) 


o 

where 

h=\i  l+^t),  s*=|(  1+rm,),  ft=|(l  +  ^t)  (A.5) 

In  this  way,  the  A-th  node  has  local  co-ordinates  rk=sk=tk=  1  in  the  rk,  sk,  tk  local 
system  of  co-ordinates. 


8 


y 

Fig.  A.l.  The  mapping  M  transforms  the  cube  into  the  hexahedron.  The  local  co-ordinates  of  nodes 
1  and  7  are  (-1  ,-l ,-l)  and  (1,1,1),  respectively. 


8 


y 


Fig.  A.2.  The  degrees  of  freedom  in  the  edge  elements  are  associated  with  the  oriented  edges  of 
the  finite  elements  mesh 
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8 


y 


Fig.  A.3.  The  (local)  numbering  scheme  for  the  faces. 

The  degrees  of  freedom  of  the  edge  elements  are  associated  with  the 
tangential  components  along  the  edges  of  the  elements  of  a  the  vector  field  to  be 
represented  (see  fig.  2)  [3]: 

U=J>,N,  (A.6) 

e=l,E 

The  shape  function  associated  with  the  edge  e={i,  j).  connecting  the  z-th  and  j'-th 
nodes  along  £,  is1  [4]: 

Ne  =^(1+Tm.)a  +  «.)V5  =seteVre  (A. 7) 

O 

where  r\e  and  t>e ,  equal  to  ±  1 ,  are  the  local  co-ordinates  of  the  edge  e  and  re,  se,  te 
is  the  local  frame  where  the  edge  e  is  described  as  se=te=  1,  re  e  [0,1] .  The  shape 
functions  for  edges  running  along  r\  and  C,  are  obtained  by  cyclic  permutations 
while  Nc,  =-Ne  when  e’={j,  i}. 

We  notice  that: 

■  The  line  integral  of  Nr  along  the  edge  e  (from  node  i  to  node  j)  is  one: 

|N„dl  =  l  (A. 8) 

<u  i 

■  The  line  integral  of  NP  is  zero  along  any  other  edge  { /,  k }  ^  {i,  j}  and 
{k,l}*{i,j}: 

1  The  ordered  pair  [i,  j }  represents  the  oriented  edge  of  the  finite  elements  mesh  directed  from 
node  i  to  node  j. 
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Appendix  B:  Tree- Cotree  Decomposition 


Solenoidality  of  the  electric  current  dipole  density  P  can  be  assured  by 
introducing  a  vector  potential  U,  defined  by 


P  =  VxU  (B.l) 

To  assure  the  unicity  of  U,  an  additional  condition  must  be  imposed.  We  prescribe 
the  component  of  the  vector  potential  along  the  direction  determined  by  an 
arbitrarily  chosen  non-vanishing  vector  field  w  which  does  not  possess  closed 
field  lines.  The  role  played  by  this  gauge  condition  is  to  eliminate  all  the 
irrotational  fields  from  the  functional  space  in  which  we  search  for  the  vector 
potentials.  In  other  words,  the  gauge  condition  should  constrain  the  null  space  of 
the  discrete  curl  operator  to  be  void  [5]. 

We  consider  a  simply  connected  region  Vd  and  we  express  U  and  P  using 
edge  elements  and  face  elements,  respectively: 

U(r)=  2>,N»  (B.2) 

e=l  ,E 

and 

P(r)=  £P,S,(r)  (B.3) 

f=lF 

where  N/s  are  edge  element  shape  functions  and  S/s  are  face  elements  shape 
functions;  consequently  Ue  and  Pf  arc  the  degrees  of  freedom  associated  with  the 
edge  e  and  the  face/. 

From  (B.l)  and  (B.2)  we  have: 

P=££/,VxN,  (B.4) 

e=l  ,E 

and  hence 

p='L'Lcf.u.sf  <b.5) 

e=l.Ef=l,F 

Equation  (B.5)  shows  that  the  coefficients  of  (B.2)  and  (B.3)  are  linearly  related: 
Cu  =p  (B.6) 

where  p  and  u  are  the  column  vectors  of  coefficients  //  and  Ue,  respectively. 

The  relationship  (B.6)  can  be  regarded  as  a  linear  system  where  u  and  p 

play  the  roles  of  unknowns  and  right  hand  side,  respectively. 

The  existence  of  a  solution  for  system  (B.6)  is  subject  to  the  compatibility 
condition,  which  requires  that  p  correspond  to  a  solenoidal  vector  field.  If  the 
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Appendix  C:  Integration  of  the  hyper- singular  term 

Here  equation  (3.15)  is  proven. 

First  we  notice  that  the  hyper- singular  term  k  ~2VVg0(r  I  r')  is  a  distribution  (see 
[9]).  This  distribution  acts  on  vector  fields  in  the  following  way: 

-^Wg0(r  I  r'):  w(r)  -»  -^v{v '  [{ So(r  1  r>(r)drj  (C.l) 

The  definition  of  the  elements  of  L„  is: 


got rlr’)-P,(r’)dr'  dr.  (C.2) 

Then,  taking  into  account  that  V  ■  P;  =0  and  applying  the  divergence  theorem 
w.r.t.  the  outer  most  integral,  it  turns  out  that: 


&0(r  I  r')-P7(r')dr'  dr 
g0(rlr’)-P;(r')dr'  >dr 
=  jj  |  n(r)  ■  P,.  (r)V  ■  J  g0  (r  I  r ')  ■  P .  (r')dr' 

*  dnf  nf 

=  TT  §  «(r)  ■  P<  (r)  I V  ■  b 0  (r  I  r')  •  p;  (r')]dr'  dS  (r ) 

^  nf 

where  n(r)  stands  for  the  outward  normal  at  point  r  of  Q.  . 

Then,  we  notice  that1 

v  Mr  I  r')P,  (r')J  =  -V'[g0(r  I  r')P,  (r')|  (C.4) 


as  follows  by  taking  into  account  the  solenoidality  of  P;  and  that  g0  is  a  function 

of  r-r’.  Finally,  replacing  (C.4)  in  (C.3)  and  applying  the  divergence  theorem 
w.r.t.  the  inner  most  integral  we  get: 


1  As  usual,  the  symbol  V  refers  to  the  spatial  co-ordinates  r  while  V’  refers  to  the  spatial  co¬ 
ordinates  r’. 
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(C.5) 


j^o 

k2 


1 

m/ 


J n(r)  ■  P,(r)  g0  (r  I  r')  n(r')  •  V,  (r')dS(r')lX(r) 

dQ.  f 
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Appendix  D:  Programming  the  finite  element  code 
List  of  symbols 

IEE  set  of  indices  of  the  cotree  edges  lying  on  the  boundary  of  element 

m 

I~E  set  of  indices  of  the  cotree  edges  lying  on  the  boundary  surface  in 

Nhf  no.  of  boundary  faces  of  the  mesh  (boundary  faces  are  numbered 

fro  1  to  Nb  j) 

Nce  no.  of  cotree  edges  of  the  mesh  (the  cotree  edges  are  numbered 

from  1  to  Nce) 

Ne  no.  of  elements  of  the  mesh  (the  elements  are  numbered  from  1  to 

Ne) 

Nedges  no.  of  edges  of  the  mesh  (the  edges  are  numbered  from  1  to  Nedges) 

Np  no.  of  nodes  of  the  mesh  (the  nodes  are  numbered  from  1  to  Np) 

Nc  no.  of  Gauss  points  per  element 

Nc  no.  of  Gauss  points  per  face 

rf'"  //-th  Gauss  point  in  the  element  m 

rnm  /7-th  Gauss  point  on  face  in 

£jTi  surface  of  the  m  -th  boundary  face  of  the  mesh 

Tm  three  dimensional  domain  occupied  by  element  m 

volfxj  volume  of  element  m 

w„  Gauss  weight  at  the  zz-th  point  of  the  master  element  [-1,  l]x[-l, 

l]x[-l,  1] 
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Procedure  2  :  compute  V 

set  V=0; 

for  m=l,  ...,  Ne 

for  n=l,  Nq 

for  ke  /f 

(m.k.n)w,  ,vo/(0; 

(y(,|i=(v(,|)t+v„,,; 

Procedure  3  :  compute  R1 2 

set  R=Q; 

for  m=l,  Ne 

for  n=l,  ...,  Nc 

for  ke  /“ 

m 

for  j  e  I(J:  and  j  <k 

P n,m,kj  =  SHPF(m,k,n)-  SHPF (m:k,n)wrvol(x iH )/ 8; 

(s Xj  =  (s)y  +  P n,m,k,j  ; 

set  (l)*  =  j  for  j  <k  ; 

set  R  = - i - R; 

c  f  —  o0 

Procedure  4  :  compute  L 

set  L=fi; 

for  m=l,  Ne 

for  1=1,  Ne 

compute  and  store  any  quantity  that  depends  only  on  m  and 

l3; 

for  ke  I 

for  je  lfE  and  j<k 

call  compute  amJ  k  . ; 

(El;  _ 

set  (L)//(  =  (L),.  for  j<k; 

Notice  that  procedure  4  is  valid  for  computing  5L  and  L  also.  In  the  first  case 
amJ  k  j  has  to  be  replaced  by  5am ,  k  .  defined  as 

1  The  matrix  R  is  sparse. 

2  As  example,  compute  and  store  G"  (r™  I  r‘.)  for  i,i’=l,  ...,NG 

3  The  algorithm  for  computing  a  m ;  k  .  is  left  to  the  reader  (see  [8]  for  the  case  of  static  scalar 
Green  function). 
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Preface 


Following  an  initial  phase  of  work  in  which  a  simulation  of  eddy  current  corrosion 
detection  for  idealized  structures  was  developed,  this  report  covers  Phase  II  of  the 
project  in  which  the  objective  is  refine  the  numerical  model  to  simulate  inspections 
taking  into  account  complex  structural  features  found  on  aircraft. 

In  the  foundation  phase  of  the  work,  the  main  objective  was  to  produce  a  mea¬ 
surement  model  for  basic  layered  structures  of  infinite  transverse  extent  containing 
local,  regular  or  irregular  flaws  representing  material  loss  due  to  corrosion.  The 
model  determines  the  magnetic  field  at  any  point  above  the  metal  structure  and  can 
therefore  predict  measurements  made  by  giant  magneto  resistive  (GMR)  sensors  or 
Hall  devices.  The  output  of  the  model  was  validated  by  comparisons  with  experi¬ 
mental  measurements  of  probe  impedance  and  sensor  signals  due  to  material  loss  in 
conducting  plates.  The  objectives  of  the  first  phase  were  achieved.  In  addition,  a 
study  was  made  of  a  specific  array  probe  concept  based  on  the  racetrack  coil.  The 
coil  shape  is  similar  in  plan  view  to  that  of  a  running  track  with  semi-circular  bends 
connected  by  straight  sections.  The  work  on  the  racetrack  coils  is  described  in  an 
article  submitted  to  Journal  of  Applied  Physics  and  is  reproduced  in  this  report. 
Also  in  the  foundation  phase,  initial  steps  were  taken  to  investigate  the  effectiveness 
of  sensor  array  probes  for  the  detection  of  roughness  on  hidden  surfaces. 

Eddy  current  measurements  on  aircraft  give  large  indications  due  to  fasteners  and 
other  structural  features  such  as  the  edge  of  a  plate  in  a  lap  joint.  It  is  often  the  case 
that  flaw  detection  is  hampered  because  the  effects  of  normal  structure  masks  the 
signal  due  to  the  defect.  Because  of  their  significance  in  the  detection  and  evaluation 
of  flaws,  the  simulation  has  been  extended  to  take  account  of  the  two  basic  structural 
features,  fasteners  and  lap  joints.  In  order  to  produce  a  reliable  simulation  code  for 
this  task,  a  new  numerical  scheme  has  been  developed  in  partnership  with  Dr’s 
Guglielnro  Rubinacci  and  Antonello  Tamburrino  from  the  University  of  Cassino,  in 
Italy. a 

This  partnership  enable  us  to  build  on  our  combined  experience  of  two  approaches 
to  the  design  of  numerical  code.  The  Cassino  group  has  been  prominent  in  the  use 
of  edge  elements  techniques  and  Dr  John  Bowler  at  CNDE  has  developed  numerical 
techniques  based  on  dyadic  Green’s  function  formulations.  Features  of  these  two 
approaches  have  been  forged  together  to  produce  a  new  formulation  with  the  aim 
exploiting  the  advantages  of  each. 

In  parallel  with  the  design  of  a  new  numerical  scheme,  experiments  have  been  car¬ 
ried  out  to  validate  model  predictions.  The  measurements  are  made  using  probes 
consisting  of  a  calibrated  Hall  sensor  located  on  the  axis  of  a  circular  drive  coil. 
Although  this  is  not  an  array  probe,  it  has  the  advantage  that  it  can  be  well  charac¬ 
terized  for  modelling  and  code  validation  purposes.  The  aim  of  the  validation  work 
is  to  check  that  the  simulation  predicts  accurate  results  for  probe  signals  due  to  an 
isolated  fastener  and  an  elementary  lap  joint. 

“Cassino  is  known  for  the  famous  World  War  II  battle  of  Monte  Cassino  which  preceded  the  fall 
of  Rome. 
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Article  1 

Simulation  of  Eddy  Current  Array  Probes:  Predic¬ 
tions  and  Measurements 

1.1  Introduction 

In  this  project  we  have  developed  a  computer  simulation  of  eddy-current  inspection  of  lay¬ 
ered  structures  using  probes  containing  giant  magneto-resistive  (GMR)  or  other  solid  state 
sensors  for  detecting  material  loss  and  surface  roughness  due  to  corrosion.  Experiments 
have  been  performed  to  validate  the  simulation  by  comparing  predictions  of  the  magnetic 
field  with  experimental  measurements  performed  using  calibrated  Hall  devices.  The  field 
calculation  treats  the  sensors  as  independent  and  therefore  only  requires  validation  using 
single  sensor  probes. 

Conventional  driver  pick-up  eddy  current  probes  use  an  induction  coil  to  induce  current 
in  the  part  and  a  pick-up  coil  to  sense  perturbation  in  the  field  due  to  flaws.  Having  a 
magnetic  field  sensor,  such  as  a  magneto-resistor  (MR)  or  Hall  device,  instead  of  the  pick¬ 
up  coil  means  that  the  sensitivity  of  the  field  measurement  is  maintained  at  low  frequency 
and  the  spatial  resolution  is  high  because  the  field  sensitive  area  in  the  sensor  is  usually 
less  than  a  millimeter  across.  There  may  be  more  intrinsic  noise  generated  by  a  solid  state 
sensor  than  by  a  pick-up  coil  but  intrinsic  noise  is  negligible  compared  with  extrinsic  noise 
due  to  surface  roughness,  material  variations,  scanning  irregularities  and  so  on.  The  main 
benefit  of  solid  state  sensors  is  the  low  frequency  sensitivity  leading  to  improved  detection 
of  subsurface  flaws. 

In  this  article,  numerical  predictions  of  magnetic  field  sensor  measurements  are  sum¬ 
marized  and  calculations  are  compared  with  experimental  measurements  to  validate  the 
simulations.  The  conclusion  reviews  the  need  to  quantify  corrosion  damage  using  array 
probe  data.  Specifically  through  assessment  of  material  loss  and  the  measurement  of  sur¬ 
face  roughness.  This  implies  that  new  inversion  tools  are  needed.  In  addition,  it  is  evident 
that  the  images  of  corrosion,  do  not  necessarily  match  up  to  the  shape  of  the  corroded  area. 
Hence  there  is  a  need  for  image  processing  tools  to  improve  the  visualization  of  damaged 
regions.  Article  2  makes  a  start  on  addressing  this  issue  by  showing  how  a  corroded  region 
may  be  imaged  quantitatively  using  traditional  signal  processing  methods.  In  this  article, 
a  brief  report  is  given  of  some  specific  numerical  results  and  validation  experiments  for  loss 
regions,  holes,  fasteners  and  lap  joints. 

1.2  Numerical  Predictions  and  Experimental  Results 

Experiments  have  been  performed  using  probes  consisting  of  a  circular  coil  and  a  Hall  sensor 
located  on  the  coil  axis,  Figure  1.  The  sensor  measures  the  magnetic  field  component  in  the 
axial  direction,  perpendicular  to  the  surface  of  the  conductor.  Two  types  of  sensor  types 
are  used  in  the  probes.  One  is  a  Honeywell  silicon  device,  634SS2,  in  a  dual  in-line  plastic 
package  which  contains  a  preamplifier.  The  other,  made  by  Asahi  Kasai  [1],  an  HW108A, 
has  a  indium  antinomide  (InSb)  sensor  [2]  which  is  more  effective  that  silicon  due  to  its 
high  carrier  mobility  but  the  package  does  not  contain  a  preamplifier.  Because  of  its  small 
size,  see  Figure  2,  the  Asahi  Kasai  seems  to  be  more  suitable  for  the  manufacture  of  an 
array  by  using  off-the-shelf  components.  Results  for  the  HW-108A  will  be  given  in  future 
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Figure  1:  Coil  and  axially  mounted  magnetic  field  sensors. 


Figure  2:  Probe  cross  section  showing  key  probe  parameters. 


reports.  Here  we  focus  on  the  measurements  taken  with  the  Honeywell  device,  which  has 
been  in  production  for  since  the  mid  eighties  and,  incidently,  is  used  in  the  probes  for  the 
Nortec  Eddyscan. 

TABLE  1.  Probe  Parameters 


Probe 

PI 

P2 

Coil  O/D  2 a  (mm) 

23.34 

28.29 

Coil  I/D  2b  (mm) 

7.86 

12.6 

Axial  length  h  (mm) 

5.0 

5.0 

Number  of  turns  N 

1760 

1750 

Coil  lift-off  zc  (mm) 

0.55 

0.55 

Resistance  (Ohms) 

128.1 

166.7 

Hall  sensor  type 

HW-108A 

634SS2 

Manufacturer 

Asahi  Kasai 

Honeywell 

Sensor  width  (mm) 

1.3 

5.3 

Sensor  length  (mm) 

2.1 

5.3 

Sensor  lift-off  zh  (mm) 

0.1 

0.1 

Hall  sensors  are  used  here  for  probe  validation  for  several  reasons,  the  two  main  ones 
being  that  the  sensor  is  linear  and  the  region  of  sensitivity  can  be  can  be  precisely  located  at 
the  Hall  slice.  In  a  GMR  for  instance,  the  sensitive  region  extends  to  the  flux  concentrator 
to  an  extent  that  is  difficult  to  quantify. 

Precise  ac  signal  measurements  are  often  plagued  by  spurious  inductive  pick-up  but  in 
the  case  of  the  Asahi  Kasai  sensor,  the  bias  current  to  the  Hall  sensor  can  be  reversed  thus 
inverting  the  phase  of  the  Hall  signal  whereas  inductive  pick-up  will  retain  the  same  phase 
following  the  bias  reversal.  By  subtracting  the  signals  measured  with  forward  and  reverse 
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Appendix  B 

Probe  Characterization  and  Calibration  Results 

Results  given  in  this  appendix  refer  to  the  characterization  and  calibration  of  probe  P3 
(see  Table  1  or  the  main  text  for  its  nominal  dimensions)  following  the  procedure  given  in 
Appendix  A. 

B1  Coil  characterization 

Having  wound  the  coil  for  probe  P3,  the  impedance  was  measured  in  air  as  a  function 
of  frequency,  Figure  Bl.  The  coil  resistance  varies  little  with  frequency,  although  at  high 
frequency  it  increases  due  to  the  skin  effect  in  the  windings. 

Typically  the  predicted  self  inductance  of  the  probe  calculated  from  the  measured 
dimensions  is  less  the  experimental  value.  This  is  attributed  to  the  fact  that  direct  mea¬ 
surement  of  the  outer  diameter  overestimates  the  effective  value  due  to  unevenness  in 
the  windings.  A  reduced  value  of  the  outer  diameter  is  adopted  for  which  the  predicted 
self-inductance  matches  the  measured  value  in  the  low  frequency  limit,  found  by  extrap¬ 
olation  of  the  measurements  to  zero  frequency.  The  low  frequency  limiting  value  of  the 
self-inductances  is  0.054634  rnH  and  the  adjusted  outer  radius  is  13.518  mm. 

Table  Bl:  Coil  Parameters 


Parameter 

Value 

Notes 

Inner  radius 

6.300  mm 

Measured 

Outer  diameter 

13.518  mm 

Fit  of  self  inductance  predictions 

Number  of  turns 

1750 

Counted 

Height 

4.860  mm 

Measured 

Liftoff 

0.710  mm 

Fit  of  half-space  impedance  predictions 

Next,  multi- frequency  impedance  data  were  obtained  for  the  probe  on  a  thick  copper 
plate.  Before  these  data  can  be  used  for  comparison  with  theory,  a  correction  must  be 
made  for  the  effects  of  parallel  capacitance  in  the  probe  cable  and  between  the  windings. 
The  Dodd  and  Deeds  [3]  model  does  not  account  for  the  effects  of  stray  capacitance  and 
so  simple  circuit  theory  is  used  in  a  compensation  scheme. 

The  experimental  impedance  change  due  to  a  thick  copper  plate  is  compared  with  Dodd 
and  Deeds  predictions  [3]  of  the  impedance  change  due  to  a  copper  half-space  conductor. 
The  coil  lift-off,  the  distance  between  the  base  of  the  coil  and  the  surface  of  the  conductor, 
is  adjusted  for  best  agreement.  Matching  with  Dodd  and  Deeds  predictions  seems  to  be  the 
most  accurate  method  of  estimating  the  coil  lift-off,  direct  measurement  being  difficult.  In 
this  case  a  value  of  0.71  mm  was  determined  by  impedance  fitting.  See  Figures  B2  for  the 
comparison  between  theory  and  experiment  for  the  coil  impedance  on  a  copper  half-space. 
The  coil  parameters  found  following  the  impedance  measurements  are  summarized  in  Table 
Bl. 


lm{Va/la}  (Amps)  Re{Va/la}  (Amps) 


e  B3:  Coil  current  Ia  and  magnetic  field  sensor  voltage  Va  variatior 
ace. 
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:  The  ratio  of  sensor  voltage  Va  to  coil  current  Ia  as  a  fun 
in  free  space. 
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Article  2 

Imaging  and  Inversion  of  Corrosion  Loss 

2.3  Estimations  of  Material  Loss 

In  eddy  current  NDE,  scan  data  is  often  presented  in  the  form  of  a  false  color  map  or 
an  ‘image’.  Typically  the  data  is  color  coded  with  minimal  pre-processing  and  because 
measurements  contain  phase  as  well  as  amplitude  information  some  of  this  information  is 
necessarily  lost  when  creating  the  image.  Ideally  one  would  like  a  map  of  the  material  loss 
rather  than  a  simple  probe  dependent  data  image  but  this  seems  to  require  a  complicated 
and  time  consuming  inversion  algorithm  which  would  extract  quantitative  material  loss 
data  from  sensor  signals  (such  a  scheme,  incidently,  has  not  yet  been  implemented  even 
though  the  methods  needed  to  accomplish  the  task  are  known  [1]). 

Image  processing  techniques  using  Fourier  transformations  and  filtering  are  fast,  simple 
and  can  reduce  the  noise  but  seem  better  suited  to  the  important  but  somewhat  cosmetic 
task  of  improving  the  appearance  of  the  image  rather  than  extracting  quantitative  estimates 
of  material  loss.  Can  we  get  a  contour  or  false  color  map  of  the  material  loss  via  a  method 
that  is  fast,  like  the  standard  image  processing  methods  and  avoids  the  computational 
burden  of  a  full  blown  inversion?  For  small  amounts  of  material  loss  in  a  flat  plate,  the 
answer  appears  to  be  “Yes” . 

2.3.1  Eddy  current  imaging 

Material  loss  in  the  form  of  an  overall  plate  thinning  can  be  estimated  by  elementary 
calibration  techniques.  One  simply  carries  out  measurements  of  the  eddy  current  signal  for 
a  number  of  different  plate  thickness  and  records  the  relationship  between  signal  and  plate 
thickness,  possibly  in  the  look-up  table.  This  is  a  very  limited  way  of  estimating  material 
loss  for  a  number  of  reasons.  For  example,  it  cannot  deal  with  local  variations  of  thickness 
and  gives  no  indication  of  the  surface  roughness.  In  addition,  it  will  not  work  were  a  signal 
due  to  a  structural  feature  interferes  with  the  material  loss  signal. 

The  last  of  these  problems  is  difficult  to  overcome  and  will  not  be  addressed  here  but 
the  other  two  can  be  dealt  with  using  a  relatively  simple  method  based  a  deconvolution 
of  the  probe  and  flaw  functions.  For  a  flat  plate,  or  for  a  stack  of  flat  plates,  with  local 
thinning  in  one  plate  due  to  say  corrosion,  the  flaw  signal  can  be  expressed  as  a  convolution 
of  a  flaw  and  a  probe  function  having  the  form 

/OO  r  OO 

/  H(x  -  x' ,y  -  y')F{x' ,y')dx'  dy  ,  (2.2) 

-OO  J  —  OO 

where  V  represents  the  measured  signal,  H  is  the  point  spread  function  of  the  probe  and 
F  the  flaw  function.  This  expression  is  valid  if  the  material  loss  is  small  compared  with 
the  plate  thickness  (say  less  than  15%).  Equation  (2.2),  may  be  written  as 

V  =  H*F  (2.3) 

where  the  *  represents  a  two  dimensional  convolution.  The  deconvolution  can  be  carried  out 
to  find  the  flaw  function  using  Fourier  transform  techniques  and  the  Weiner  filter.  Thus, 
an  approximation  of  the  deconvolved  function  is  obtained  from  the  inverse  2-D  Fourier 
transform  of 
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—  J-f*  — 

F= ^ - V  (2.4) 

HH*  +  I< 

where  F  is  the  2-D  Fourier  transform  of  the  flaw  function,  F  is  the  2-D  Fourier  transform  of 
the  point  spread  function  and  K  is  a  constant  depending  on  the  noise  in  the  data.  Equation 

(2.2)  is  justified  later  in  this  article  but  first,  it  is  helpful  to  recall  its  previous  use  in  eddy 
current  signal  processing. 

2.3.2  Surface  crack  imaging 

The  idea  of  expressing  an  eddy  current  probe  signal  as  a  convolution  is  not  new.  An  early 
attempt  to  deconvolve  eddy  current  signals  due  to  surface  defects  such  as  cracks  in  turbine 
disks  was  made  by  Joyson,  McCary,  Oliver  Silverstein-Hedengren  and  Thumhart  [2]  at 
GE’s  Corporate  Research  and  Development.  Clearly,  the  work  at  GE  and  elsewhere  which 
treats  eddy  current  images  as  convolutions,  would  have  a  better  foundation  if  it  could  be 
justified  theoretically  for  surface  cracks.  Equation  (2.2)  is  simply  assumed,  as  a  working 
hypothesis,  in  the  hope  that  the  results  justify  the  approach  d  posterior  but  the  physics  of 
the  problem  seems  at  variance  with  the  assumption. 

The  interaction  with  a  crack  depends  on  the  direction  of  the  induced  current  which 
must  be  taken  into  account  even  if  the  formula  used  for  the  measured  signal,  equation 

(2.2) ,  is  a  scalar  formula.  In  the  GE  work,  a  tape  head  probe  was  used  to  control  the 
field  direction  giving  a  flaw  image  which  is  dependent  on  the  probe  orientation.  There  are 
disadvantages  in  being  restricted  to  tape  head  probes  and  having  to  deal  with  crack  images 
that  are  orientation  dependent.  However,  the  probe  response  to  subsurface  corrosion  is 
not  strongly  direction  dependent,  which  suggests  that  deconvolution  in  this  context  may 
successful  for  a  wide  range  of  probes. 

A  tinre-harnronic  field  is  typically  represented  by  a  phasor,  hence  the  interaction  with 
a  surface  crack  is  described  by  a  complex  vector  process  in  three  dimensions.  It  is  not 
clear  how  the  crack  signal  would  reduce  theoretically  to  a  real  scalar  convolution  in  two 
dimensions  as  represented  by  (2.2).  However,  the  theoretical  reduction  to  a  complex  2D 
scalar  convolution  can  be  performed  for  limited  subsurface  corrosion  loss  in  a  plate. 

2.3.3  Imaging  subsurface  corrosion  in  a  plate 

Although  (2.2)  has  not  been  justified  theoretically  for  surface  flaws  such  as  cracks,  it  can 
be  supported  for  small  amounts  of  material  loss  in  plates.  The  reason  for  this  is  that  plate 
thinning  gives  a  small  perturbation  of  the  field  over  a  relatively  large  region  whereas  a 
crack  gives  a  large  localized  field  perturbation. 

In  general  the  flaw  signal  detected  by  a  driver  pick-up  probe  can  be  derived  rigorously 
from  a  reciprocity  principle  [3]  to  give,  suppressing  the  probe  coordinates, 

V] 2  =  -  /  E(1)(r')  •  P(r')  dr'  (2.5) 

J  n0 

for  unit  coil  current,  where  Do  is  the  flaw  region,  E^^r')  is  the  unperturb  field  at  the 
flaw  due  to  the  pick-up  sensor2  and  P(r)  is  the  current  dipole  density  of  the  flaw  induced 

2If  the  pick-up  happens  to  be  a  magnetic  field  sensor  then  the  corresponding  electric  field  is  that  produced 
by  a  magnetic  dipole  at  the  site  of  the  sensor. 
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and  summing  to  obtain  the  total  field.  Section  II  gives  the  theory  for  the  straight  linear  coil 
elements  and  Section  III  describes  the  solution  for  the  D-coil  representing  a  semicircular 
bend.  The  results  and  conclusions  follow  the  analysis  sections. 

3.2  Linear  Coil  Field 

The  following  account  describes  the  calculation  of  the  electric  field  induced  in  a  conductor 
by  a  tinre-harnronic  current  in  a  linear  coil  consisting  of  only  the  straight  parts  of  the 
racetrack  coil  shown  in  Fig.  1.  The  current  path  is  closed  by  joining  the  ends  of  the 
straights  by  current  filaments  but  ultimately  the  effect  of  these  filaments  is  cancelled  by 
similar  filaments  added  to  the  representation  of  the  bends.  Results  are  given  for  a  field  in 
a  homogeneous  half-space  conductor  in  the  region  z  <  0,  and  an  infinite  conducting  plate. 
However,  similar  results  for  layered  conductors  are  readily  obtained  by  simply  changing  the 
Green’s  functions  used  in  the  present  calculation  in  favor  of  one  which  embodies  the  correct 
interface  conditions  of  a  stratified  conductor.  The  magnetic  dipole  formulation  used  in  this 
section  to  represent  the  field  of  the  straights  is  also  used  for  the  bends  in  the  next  section. 

Consider  a  non-magnetic  conductor  occupying  the  half-space  defined  by  z  <  0,  excited 
by  a  current  source.  The  electric  field  in  adjoining  half-spaces  is  a  solenoidal  solution  of, 

V2E(r)  =  jw/io  J  (r),  z  >  0  and  (V2  —  ju/iocrj  E  (r)  =  0,  z  <  0,  (3.1) 

where  a  is  the  conductivity  of  the  conductor.  The  electric  field,  being  transverse  to  the 
z— direction  and  having  zero  divergence,  can  be  expressed  in  terms  of  a  transverse  electric 
(TE),  scalar  potential: 

E  (r)  =  -jujfj, oV  x  zip'( r),  (3.2) 

where  z  is  a  unit  vector  in  the  preferred  direction.  The  transverse  source  current  J(r), 
having  zero  divergence,  can  similarly  be  written  in  transverse  scalar  form  as 

J(r)  =  —V  x  [zM(r)j.  (3.3) 

V  o 

The  function  M(r)  represents  the  current  source  in  terms  of  the  magnetic  dipole  density, 
the  orientation  of  the  polarization  being  in  the  z— direction.  This  is  an  adaptation  of  the 
magnetic  shell  model  which  represents  a  filamentary  current  loop  in  terms  of  magnetic  shell 
bounded  by  the  loop.  Here,  the  magnetic  dipole  distribution  occupies  a  volumetric  region 
between  the  upper  and  lower  extent  of  the  coil  where  h  +  c  >  z  >  h  —  c,  2c  being  the  height 
of  the  coil  and  h  the  height  of  the  mid  point  of  the  coil  above  the  surface  of  the  conductor. 
Equations  (3.2)  and  (3.3)  are  substituted  into  (3.1)  to  give 

V2V/(r)  =  —  — -M(r),  z>0  and  (v2  —  y/(r)  =  0,  z  <  0.  (3.4) 

Vo  v  7 

An  expression  for  the  solution  in  terms  of  a  Green’s  function,  satisfying 

V2G(r,  r')  =  —  e)(r  —  r'),  z>0  and  (V2  —  jui/j.Qo'j  G(r,  r7)  =  0,  z  <  0,  (3.5) 


i  r 

V>'(r)  =  - 


is  written  as 


1 


G( r,  rr)  M(r')  dr, 


(3.6) 
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3.4  Results 

The  sensor  signals  due  to  a  square  recess  in  the  bottom  surface  of  a  plate  of  thickness  4.85 
mm,  Fig.  2,  with  the  racetrack  coil  above  the  plate  providing  the  excitation  field  has  been 
calculated  using  a  volume  element  code  [4].  The  dimensions  of  the  recess  are  25.4  mm  x 

25.4  mm  x  3  mm.  These  and  other  parameters  are  given  in  Table  1. 

The  sensors  measure  the  magnetic  field  component  normal  to  the  surface  of  the  con¬ 
ductor.  Fig.  3,  compares  the  variation  with  probe  position  of  the  normal  magnetic  field 
at  a  central  sensor.  Field  values  are  normalized  to  a  coil  current  of  1  Amp.  The  field 
variation  is  due  to  the  back  surface  recess  is  plotted  for  three  racetrack  excitation  coils. 
One  has  32  mm  straight  sections,  one  has  16  mm  straight  sections  and  the  other  has  no 
straight  sections  and  is  thus  a  circular  coil.  Results  from  the  zero  straight  section  coil 
using  computer  code  for  the  racetrack  analysis  agree  with  results  for  a  dedicated  circular 
coil  calculation  [5] .  The  absolute  value  of  the  z-component  of  the  field  found  by  simulating 
the  response  of  a  33  element  sensor  array  is  shown  in  Fig.  4. 


3.5  Conclusion 

The  theory  for  a  racetrack  coil  in  the  presence  of  a  stratified  conductor  has  been  given 
in  two  parts  based  on  a  formulation  using  a  magnetic  dipole  representation  of  the  effect 
of  the  coil.  In  the  first  part  the  field  due  to  the  straight  sections  of  the  track  are  found 
using  a  two  dimensional  Fourier  transform,  and  in  the  second  part,  the  field  due  to  the 
semicircular  bends  of  the  track  are  determined  using  integrals  containing  Bessel  functions. 
The  racetrack  coil  geometry  will  be  incorporated  into  a  probe  design  in  which  an  array 
of  magnetic  sensors  are  located  along  the  center  line  of  the  track  parallel  to  the  straight 
sections.  In  this  way  the  local  applied  field  experienced  by  each  sensor  is  similar  yet  the 
probe  itself  is  compact  and  easy  to  manipulate. 
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FIG.  3.  Variation  with  probe  position  of  the  magnitude  and  phase  of  the  magnetic  field  at 
the  site  of  a  central  sensor  0.89  mm  from  the  surface  of  the  material. 
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Article  4 

Dyadic  Green’s  Functions  for  the  Calculation  of 
Eddy  Currents  in  Planar  Stratified  Conductors 

4.1  Introduction 

The  use  of  dyadic  Green’s  functions  for  electromagnetic  field  problems  is  discussed  in  refer¬ 
ence  [1].  The  dyadic  Green’s  function  for  an  infinite  stratified  conductor  in  the  quasi-static 
limit  can  be  derived  using  scalar  decomposition.  The  approach  [2]  reduces  the  problem  to 
one  of  finding  scalar  transverse  magnetic  (TM)  and  transverse  electric  (TE)  components 
of  the  field  in  a  multi-layered  planar  structure.  Here  we  give  specific  results  for  a  half¬ 
space  conductor  and  an  infinite  slab.  The  tedious  but  straightforward  steps  describing  the 
derivation  of  the  TE  and  TM  scalar  Green’s  functions  satisfying  the  appropriate  continuity 
conditions  at  layer  interfaces  are  omitted. 


4.2  Dyadic  Green’s  Function 


4.2.1  Scalar  Decomposition 


The  required  dyadic  Green’s  function,  representing  physically  the  electric  field  due  to  a 
singular  electric  source,  can  be  written  in  the  form  [2]  [3] 

Q{r\r')  =  p(V  x  V  x  z)(V  x  V'  x  z)U'( r|r')  +  (V  x  z)(V'  x  z)U"( r|r').  (4.1) 

where  the  ‘wavenumber’  k  having  a  positive  real  part  is  found  from  k 2  =  icofra.  This 
scalar  form  is  valid  everywhere  except  at  the  singularity  [2]  but  this  need  not  concern  us 
because  the  correct  representation  of  the  singularity  is  well  know  and  in  any  case  will  be 
reintroduced  later  [by  using  eq  (4.8)].  The  numerical  treatment  of  the  hyper-singularity 
will  not  be  discussed  here  however  reference  [4]  covers  this  topic  adequately. 

The  functions  U'  and  U"  are  related  to  transverse  magnetic  (TM)  and  transverse  electric 
(TE)  Green’s  functions  by 


G"(r|r') 

G"(r|r') 


=  -v? 


U'  (r|rO 
U"(  r|r') 


(4.2) 


In  homogeneous  regions,  the  scalar  Green’s  functions,  G'  and  G" ,  satisfy  the  inhomogeneous 
Helmholtz  equation: 


(V2  +  k2) 


G\r\r') 
G"(r  |r') 


5(r  -  r') 


(4.3) 


subject  to  appropriate  interface  conditions  which  are  derived  from  the  continuity  of  the 
tangential  electric  and  magnetic  field. 


4.2.2  Interface  Conditions 

Specifically,  the  continuity  of  the  tangential  electric  field  means  that3 

3The  square  bracket  is  used  to  denote  the  jump  of  a  function  at  an  interface.  Thus  the  expression  [/]  =  j 
means  that  the  jump  (discontinuous  change)  in  /  at  an  interface  is  j.  If  it  happens  that  [f]=0  at  an  interface 
then  it  means  that  /  is  continuous  there. 
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4.4.2  Dyadic  Forms 

To  assemble  the  explicit  slab  dyadic  Green’s  function,  we  sub  the  explicit  form  of  the  scalar 
Green’s  functions  into  equation  (4.1).  The  result  is  then  manipulated  using  (4.8)  and  (4.9) 
(see  Felsen  and  Marcuvitz  [2],  pl8)  to  give 

Q(r\r')  =  Qt[y\y')  +  £A(r|r'),  (4.22) 

where 

^x(r|r')  =  t/(°)(r|r')  +  <j(0)(r|r+_)  +  £(0)(r|r_+) 

+  W_(r|r')  +  V  x  z[V  x  zF_(r|r')]  (4.23) 

and 

Q\{y\y')  =  <?W(r|r+)  +  £w(r|r_) 

+  W+(r|r')  +  V  x  z[V'  x  zF+(r|r')].  (4.24) 

The  new  variables  are  defined  as  follows: 

r+  =  r  —  2(z'  —  d)z 
r_  =  r  -  2 (z')z 
r+_  =  r'  —  2  dz 

r_+  =  r'  +  2  dz.  (4.25) 

f?o(r|r)  is  the  dyadic  Green’s  function  for  an  unbounded  conductor  given  by 

g(0)(r|r')  =  [X+  -LVV]Go(r|r'),  (4.26) 

with 


Go(r|r')  = 


Dife|r— r'| 


(4.27) 


47r| r  —  r'|  ’ 

The  remaining  terms  in  (4.23)  are  due  to  reflection  at  the  upper  and  lower  surfaces  of  the 
source  stratum.  The  first  two  of  these  remaining  terms  are  double  refection  image  terms. 
In  (4.24)  we  find  that  two  single  reflection  image  terms  given  by 


1 


<?W(r|r±)  =  [!'  - -VV'Mrlri), 
where  Tl  =  xx  +  yy  —  zz.  The  other  terms  are  defined  as  follows: 

W_(r|r')  =  [X  +  pVV]W_(r|r'), 


W+(r|r')  =  [X,-^VV/]IT+(r|r/), 


(4.28) 

(4.29) 

(4.30) 


where 


1  fOO  ir 

W-(p,  Z,  z')  =  —  /  —  T'(«)  e7(*-*'-2d)  +  e-7(*-*'+2d)  J^Kp)KdK  (4.31) 
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Article  5 

Edge  Elements  Based  Numerical  Computation  of 
Dipole  Density  in  Eddy  Current  Testing  Problems 

Guglielmo  Rubinacci  and  Antonello  Tamburrino 

DAEIMI 

Laboratory  of  Computational  Electromagnetics 
and  Electromagnetic  Nondestructive  Evaluation 
Universita  degli  Studi  Cassino,  Via  G.  Di  Biasio,  43  Cassino,  1-03043 


1.  Introduction 


This  report  presents  an  edge  elements  based  numerical  method  for  computing  the 
dipole  density  P  arising  from  the  interaction  of  eddy  currents  with  a  flaw  in  a 
conductive  slab. 

The  dipole  density  appears  in  the  region  Q.  f  occupied  by  the  flaw  and  it  can 
be  determined  by  solving  an  integral  equation  in  Clf .  This  is  a  great  advantage 
from  the  numerical  viewpoint  because  it  requires  the  discretization  only  of  Clf , 

which  is  generally  a  small  part  of  the  whole  conductive  domain  £l0 .  However,  the 
electric-electric  dyadic  Green’s  function  in  the  presence  of  the  conductive 
material  (without  flaw)  is  required  and  the  kernel  of  the  resulting  integral  equation 
is  hyper-singular.  Moreover,  the  electric-electric  dyadic  Green’s  function  is 
known  only  for  canonical  geometries  such  as,  for  instance,  the  homogenous  slab. 

Assuming  that  the  conductivity  of  the  flaw  is  constant,  the  dipole  density  is 
solenoidal  and  can  thus  be  represented  as  the  curl  of  a  vector  potential  U,  which 
we  term  the  dipole  vector  potential.  The  numerical  formulation  is  based  on  the 
Galerkin  method  applied  to  the  integral  equation  satisfied  by  P  and  expanding  U 
as  the  linear  combination  of  edge  element  based  shape  functions.  The  uniqueness 
of  the  dipole  vector  potential  is  numerically  imposed  by  means  of  the  tree-cotree 
decomposition  of  the  finite  elements  mesh  used  to  build  the  edge  element  shape 
functions.  This  numerical  formulation  automatically  takes  into  account  the 
solenoidality  of  P  and  thus  allows  the  number  of  unknowns  to  be  reduced  with 
respect  to  the  method  of  moment  when  applied  on  the  same  mesh,  and  makes  it 
possible  to  improve  the  rate  of  convergence  of  the  numerical  solution.  Moreover, 
by  taking  into  account  the  solenoidality  of  P,  the  volume  integral  arising  from  the 
hyper-singular  part  of  the  dyadic  Green’s  function  is  reduced  to  a  surface  integral 
containing  the  scalar  Green’s  function  as  its  kernel. 
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The  first  term  in  (3.16)  corresponds  to  a  non-singular  kernel  and  its  integration 
can  be  carried  out  by  standard  numerical  integration  techniques,  whereas  the 
second  term  correspond  to  the  static  scalar  Green’s  function  and  its  integration 
can  be  carried  out  as  showed  in  [8]. 
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(A.9) 


|N,  dl  =  0 

{Wl 

■  The  tangential  components  of  N,,  are  continuous  across  the  face  of 
adjacent  elements  since  the  scalar  functions  Nk  are  continuous. 

■  The  normal  components  of  Ne  are  not  necessarily  continuous. 

■  N,,  and  VNk  belong  to  the  same  functional  space:  the  gradient  of  a  nodal 
shape  function  is  given  by  a  linear  combination  of  the  edge  shape 
functions  having  in  common  that  node: 


W,  =  £g„n, 

e=l  ,E 


(A.  10) 


where  E  is  the  number  of  edges  of  the  mesh  and  CL  is  the  ExNP  incidence 
matrix  defined  by: 


G 


em 


+  1  if  e  =  {i,  j} and 
-lif  e  =  {i,j}md 


m  =  j 
m  =  i 


(A.ll) 


[0  otherwise 


Similarly,  a  set  of  shape  functions  S/  associated  with  the  faces  of  the  elements  can 
be  introduced.  In  this  case,  the  flux  of  Sy  is  one  across  face /and  zero  across  any 
other  faces.  The  normal  component  of  S/is  continuous  across  the  faces  of  adjacent 
elements,  whereas  the  tangential  component  is  not  necessarily  continuous. 
Moreover,  it  can  be  verified  that  the  curl  of  an  edge  shape  function  is  a  linear 
combination  of  face  shape  functions: 


VxN,=  £c„  S, 

f=hF 


(A.  12) 


where  F  is  the  number  of  faces  of  the  mesh  and  £  is  the  F/E  incidence  matrix 
defined  by1: 


C 


fe 


+  lif  /  =  {m,  n,l, /} and  e  =  {m,  n} 
■  -  lif  /  =  {m,n,l,i} and  e  =  {n, m} 
0  otherwise 


(A.  13) 


Notice  that  from  (A. 8)  and  (A.9)  it  follows  that  the  edge  element  shape  functions 
N/s  are  linearly  independent.  Similarly,  the  face  element  shape  function  S/s  are 
linearly  independent. 

Moreover,  the  matrix  £L  represents  the  discretization  of  the  gradient  operator 
whereas  the  matrix  C  the  discretization  of  the  curl  operator: 

cp(r)=  £0iJV,(r)^V(p=  £(Ga)tN,  (A.14) 

k-l,NP  e=l  ,E 

U(r)=  £,N,W^VxD=  £  (Cb)  S,  (A.  15) 

e=l  ,E  f=\,F 


1  The  ordered  sequence  {m,n,l,i }  represents  the  face  having  m,n,l,i  as  nodes.  Notice  that  any 
circular  permutation  represents  the  same  face. 
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where  a  and  b  stand  for  the  column  vectors  of  numerical  coefficients  ak  and  bk. 
The  curl  of  N,,  has  the  following  expression: 

VxN  e=V(sje)xVre  (A.16) 


The  explicit  expression  for  the  gradient  is  obtained  by  the  usual  rules  of  partial 
differentiation,  using  the  Jacobian  matrix  J  of  the  transformation  M: 


Vr  = 


dr 

dr 

dr 

T 

_i 

dr 

dr. 

dr 

e 

e 

e 

=  J 

e 

e 

e 

dx 

dy 

dz 

dre 

dse 

d^_ 

where 


dx 

dy 

dz 

dre 

dr. 

dre 

dx 

dy 

dz 

dse 

dse 

dse 

dx 

dy 

dz 

1*. 

dte 

K. 

=  _ 

2xv 

k  v2  v3] 


3  V3XVl  V  !  X  V  2  ] 


=  j  '[l  0  of 


(A.  17) 


(A.  18) 


(A.  19) 


and  the  determinant  of  the  Jacobian  matrix  is: 

det(j)=  v2  x  v3  •  v,  (A. 20) 

Cyclic  permutations  give  the  other  gradients. 
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for  T  =  1,  ....  Nbf 

compute  and  store  any  quantity  that  depends  only  on  in 
and  l  1; 

for  k  e  I~E 

for  je  I~E  and  j<k 

call  compute  fL  T  k  ; 

y,=y,+iw 

We  notice  that  the  volume  and  surface  integrals  are  carried  out  on  the  master 
element  [-1,  l]x[-l,  l]x[-l,  1]  or  the  master  face  [-1,  l]x[-l,  1],  respectively.  The 
Gauss  weight  refer  to  the  master  element  or  face;  the  Gauss  points  for  a  generic 
element  or  face  are  obtained  as  image  under  the  transformation  M  of  the  Gauss 
points  for  the  master  element  and  face,  respectively. 


1  As  example,  compute  and  store  g(i (r;m  I  rj )  for  i,i’=l,  ...,  Na 
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Procedure  2  :  compute  the  incidence  matrix  ( Nedges ,  inc) 


lat=0; 

for  id=l,  ...,  Ne 

for  each  edge  i  (i=l,  ...,  12)  of  element  iet 
set  ind=false; 
if  edge  i  is  already  in  inc1 

then  set  latel(i,  iei)=position  where  edge  i  occurs  in  inc; 
else 

lat=lat+l; 

insert  the  label  of  nodes  corresponding  to  edge  i  in  inc(:, 

lat); 

latel(i,  iei)=lat; 

hi ;dges — lat , 

Procedure  3  :  compute  ltact2 
for  each  element  iel 

for  each  edge  i  (i=l,  ...,  12)  of  element  ie/ 
lat=latel(i,ie/); 

if  lat  belong  to  the  co-tree  then  set  ltact(i,  iei)=lat; 
else  ltact(i,  iei)=0; 


Procedure  4  :  compute  ltact_bou 

for  each  boundary  face  jf 

for  each  edge  i  (i=l,  ...,  4)  of  the  face  jf 
ie,=ielface(jf); 
ilat=indface(  i,jf); 
lat=ltact(ilat,iei); 
ltact _bou(  i,jf)=lat; 


1  This  test  consists  in  searching  the  column  of  inc  containing  the  nodes  related  to  edge  i. 

2  Notice  that  the  procedure  for  the  tree -cotree  decomposition  is  described  in  [6].  The  input  of  this 
procedure  is  the  incidence  matrix  inc  and  the  output  is  the  list  of  co-tree  edges  ncoalb. 
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